heaviside function in fun for fminbnd
4 views (last 30 days)
Show older comments
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{:});
0 Comments
Answers (2)
Walter Roberson
on 3 Feb 2017
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.
0 Comments
See Also
Categories
Find more on Assumptions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!