heaviside function in fun for fminbnd

I have a function

fun = @(x) f(x)*heaviside(g(x))

so it is a function of x but within it also a heaviside function that depends on x as well.

I want to find a minimum value for x by using fminbnd:

x2= fminbnd(@(x) fun(x),0,minel);

but I get a lot of errors:

Error using symengine (line 59)
An arithmetical expression is expected.
Error in sym/privUnaryOp (line 909)
            Csym = mupadmex(op,args{1}.s,varargin{:});
Error in sym/heaviside (line 14)
Y = privUnaryOp(X, 'symobj::map', 'heaviside');
Error in
@(x)abs(mtotalh2(m)-(x*pin(i)*10^6/(heaviside(minverhI-(capel(e)/((x)*pin(i))))*(c2*log10(opslagcapaciteit/(capel(e)/((x)*pin(i))))+c1)+heaviside((capel(e)/((x)*pin(i)))-minverhI)*(c2*log10(opslagcapaciteit/minverhI)+c1))/(96500*1007.5)+(1-x)*pin(i)*10^6/(heaviside(minverhI-(capel(b)/((1-x)*pin(i))))*(c2*log10(opslagcapaciteit/(capel(b)/((1-x)*pin(i))))+c1)+heaviside((capel(b)/((1-x)*pin(i)))-minverhI)*(c2*log10(opslagcapaciteit/minverhI)+c1))/(96500*1007.5)*IelI(i,e,b,m,p)))
Error in @(x)fun(x)
Error in fminbnd (line 292)
    fu = funfcn(x,varargin{:});

Answer 1

[Torsten]

Use

fun = @(x) f(x)*(g(x)>0)

Best wishes

Torsten.

[Maurice Hendriks]

Thanks for your reply! I tried this but it still gives the same errors. The only thing that makes it different from a 'normal' function is that there is a log10(h(x)).

 fun = @x f(x)*(g(x)>0)*log10(h(x))

I think the fminbnd function can't work with the log10 part because before it worked with similar functions without the log10. Any ideas how to solve this?

Answer 2

[Walter Roberson]

Be careful: the value of heaviside(0) is not formally defined as either 0 or 1. You can change the result in the Symbolic Toolbox by using symprefs(); by default it is 1/2, which is typically not what people expect.

If your h(x) can be 0 or infinity, then log10(h(x)) would be +/- infinity there. If your g(x) happens to be <= 0 at the location that h(x) is 0 or infinity, then that would give rise to 0 * infinity, which is defined as NaN

There is no numeric solution in MATLAB: you need to use a control structure.

You appear to be using the symbolic toolbox at some point. If you are doing that at for more reasons than using heaviside, then I suggest you build your formula to be optimized using piecewise() (needs R2016b or later) and then use the trick I recently found: use matlabFunction() and tell it to write the result to a file. The generated file will use "if" as appropriate. And if at all possible, do that building of the formula outside of the function you give as the objective function; for example use it to build the function handle that you pass in as your objective function.

Note: using matlabFunction() to generate "if" will not work if you are trying to vectorize, since it does not use logical indexing to implement the tests: it assumes it is working with a scalar. Also, there are bugs in the code it generates if the piecewise() expression is an element of a vector or array. Fortunately, these conditions are not a problem for fmincon.

Another note about using matlabFunction to build for fmincon: if you have multiple variables that you want to correspond to the different elements of the "x" that gets passed to the objective function, then be sure to use the 'vars' parameter of matlabFunction to pass a column vector of variable names:

fun = matlabFunction( Expression, 'vars', {[alpha; beta]}, 'file', 'objfun.m')

would be for expecting x(1) to match to alpha and x(2) to match to beta in the expression.

A bunch of what is in your formula

abs(mtotalh2(m)-(x*pin(i)*10^6/(heaviside(minverhI-(capel(e)/((x)*pin(i))))*(c2*log10(opslagcapaciteit/(capel(e)/((x)*pin(i))))+c1)+heaviside((capel(e)/((x)*pin(i)))-minverhI)*(c2*log10(opslagcapaciteit/minverhI)+c1))/(96500*1007.5)+(1-x)*pin(i)*10^6/(heaviside(minverhI-(capel(b)/((1-x)*pin(i))))*(c2*log10(opslagcapaciteit/(capel(b)/((1-x)*pin(i))))+c1)+heaviside((capel(b)/((1-x)*pin(i)))-minverhI)*(c2*log10(opslagcapaciteit/minverhI)+c1))/(96500*1007.5)*IelI(i,e,b,m,p)))

is constant in x, and should be calculated outside for efficiency. You should be left with no indexing, reduced function calls, and with common expressions that are independent of x moved out; leaving only expressions in which x is a necessary component of the calculation.