how can i avoid Nan in matlab expression and return 0
8 views (last 30 days)
Show older comments
The below three expression gives three different values,
and can someone advice me how can i avoid Nan values for 0 input and return 0 value for it
(sqrt((-220).^2)./-220+1)/2.*(10-5) + 5
(sqrt((220).^2)./220+1)/2.*(10-5) + 5
(sqrt((0).^2)./0+1)/2.*(10-5) + 5
>> (sqrt((-220).^2)./-220+1)/2.*(10-5) + 5
(sqrt((220).^2)./220+1)/2.*(10-5) + 5
(sqrt((0).^2)./0+1)/2.*(10-5) + 5
ans =
5
ans =
10
ans =
NaN
4 Comments
Walter Roberson
on 5 May 2021
NN is trying to use Problem Based Optimization for something it was never designed to do, so the straight-forward methods are not available, and instead it is necessary to find some mathematical "cheat" that the Optimization Toolbox will permit.
Answers (5)
Matt J
on 5 May 2021
Edited: Matt J
on 5 May 2021
Replace
sqrt(x.^2)./x
with
sign(x)
5 Comments
Walter Roberson
on 5 May 2021
Optimization variables cannot do logical multiplication. They can only use comparisons in the context of constraints.
Image Analyst
on 5 May 2021
Why can't you just simply assign it to a variable and check if that variable is nan, and if it is, assign it to zero?
result = (sqrt((0).^2)./0+1)/2.*(10-5) + 5;
if isnan(result)
result = 0;
end
1 Comment
Walter Roberson
on 6 May 2021
Because the context is Problem Based Optimization, which does not support if or isnan() and only permits comparisons as part of constraints.
Walter Roberson
on 6 May 2021
The following is not exactly right:
M = optimvar('M', 1, 5)
costa = 5
costb = 10
delta = eps(realmin)
part1a = (M - sqrt(M.^2))/2
part1b = part1a./(part1a - delta)
part2a = (M + sqrt(M.^2))/2
part2b = part2a./(part2a + delta)
cost = part1b * costa + part2b * costb
Mathematically it is wrong at exactly two points, M = -eps(realmin) and M = +eps(realmin) . For those two points, the output should be costa and costb respectively, but instead the formula mathematically gives costa/2 and costb/2 at those two points instead.
In practice, though, for values sufficiently close to +/- realmin, numeric evaluation might return costa+costb and for values sufficiently clost to eps(realmin) numeric evaluation might return 0 instead of costa or costb .
The exact result in a range close to +/- realmin is going to depend on the exact order of evaluation, which is not something that we have control over; optimization could potentially re-arrange the evaluation.
This code has been constructed so that it should never return NaN.
How important is it for your purposes that the values must be correct near +/- realmin, given that it is designed to return 0 for 0 exactly?
4 Comments
Walter Roberson
on 6 May 2021
Example. This needs to be upgraded to have Costb passed in as well, but your scripts and the Simulink model currently only pass in a single cost vector.
function [Pgrid,Pbatt,Ebatt] = battSolarOptimize(N,dt,Ppv,Pload,Einit,Cost,FinalWeight,batteryMinMax)
% Minimize the cost of power from the grid while meeting load with power
% from PV, battery and grid
prob = optimproblem;
% Decision variables
PgridV = optimvar('PgridV',N);
PbattV = optimvar('PbattV',N,'LowerBound',batteryMinMax.Pmin,'UpperBound',batteryMinMax.Pmax);
EbattV = optimvar('EbattV',N,'LowerBound',batteryMinMax.Emin,'UpperBound',batteryMinMax.Emax);
% Minimize cost of electricity from the grid
prob.ObjectiveSense = 'minimize';
%{
prob.Objective = dt*Cost'*PgridV - FinalWeight*EbattV(N);
%}
Costa = Cost(:,1);
Costb = Cost(:,end) + randn(size(Costa))/20;
r = 200;
%Cost = (PbattV<=r).*Costa+ (~PbattV>=r).*Costb;
Cost = fcn2optimexpr(@(PV) (PV < r) .* Costa + (PV > r) .* Costb, PbattV);
%P1 = dt*Cost'*PbattV;
%P2 = dt*Cost'*PbattV;
P = dt*Cost'*PbattV;
prob.Objective = dt*Cost'*PgridV - FinalWeight*EbattV(N);
% Power input/output to battery
prob.Constraints.energyBalance = optimconstr(N);
prob.Constraints.energyBalance(1) = EbattV(1) == Einit;
prob.Constraints.energyBalance(2:N) = EbattV(2:N) == EbattV(1:N-1) - PbattV(1:N-1)*dt;
% Satisfy power load with power from PV, grid and battery
prob.Constraints.loadBalance = Ppv + PgridV + PbattV == Pload;
x0.PgridV = ones(1,N);
x0.PbattV = batteryMinMax.Pmin * ones(1,N);
x0.EbattV = batteryMinMax.Emin * ones(1,N);
% Solve the linear program
options = optimoptions(prob.optimoptions, 'Display', 'iter', 'MaxFunctionEvaluations', 1e7, 'MaxIterations', 1024 );
[values,~,exitflag] = solve(prob, x0, 'Options', options);
% Parse optmization results
if exitflag <= 0
fprintf('warning: ran into iteration or function evaluation limit!\n');
end
Pgrid = values.PgridV;
Pbatt = values.PbattV;
Ebatt = values.EbattV;
Matt J
on 6 May 2021
Edited: Matt J
on 6 May 2021
Here's a way you can do it by adding some additional binary variables and linear constraints. It requires that x be bounded to the interval [-1,1]. You can easily introduce a variable z=A*x+B if you need a variable that spans a different range.
x=optimvar('x','LowerBound',-1,'UpperBound',+1);
b1=optimvar('b1','LowerBound',0,'UpperBound',1,'type','integer'); %Binary variables
b2=optimvar('b2','LowerBound',0,'UpperBound',1,'type','integer');
%%EDITED
con1(1)= b1>=x+eps(0); % b1 "ON" when x>=0
con2(1)= b1<=x+1; % b1 "OFF when x<0
con1(2)= b2>=x; % b2 "ON" when x>0
con2(2)= b2<=x+1-eps(1) % b2 "OFF when x<=0
y = 5*(1 + (b1+b2)/2 -3*(b1-b2)/2);
% x>0 ==> b1=b2=1 ==> y=10
% x<0 ==> b1=b2=0 ==> y=5
% x=0 ==> b1=1, b2=0 ==> y=0
3 Comments
Walter Roberson
on 7 May 2021
Each constraint is structed by constraint type (same for all elements in the constraint), and two pieces of data indicating what the constraint is operating on.
I do not know why they choose to implement that way..
See Also
Categories
Find more on Simulink Real-Time 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!