MATLAB Answers

how can i avoid Nan in matlab expression and return 0

7 views (last 30 days)
NN
NN on 5 May 2021
Commented: Walter Roberson on 7 May 2021
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
NN
NN on 5 May 2021
:-(
Yes, but how logical comparison can be performed inside optimisation problem ? I am worried if there is no solution for this

Sign in to comment.

Answers (5)

Matt J
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
Walter Roberson on 5 May 2021
Optimization variables cannot do logical multiplication. They can only use comparisons in the context of constraints.

Sign in to comment.


Stephen Cobeldick
Stephen Cobeldick on 5 May 2021
Edited: Stephen Cobeldick on 5 May 2021
"if X is positive it must give 10, If X is negative it must give 5, If X is zero, it must give 0."
X = randi([-3,3],1,9) % random data
X = 1×9
0 -3 0 1 2 -2 3 -3 0
V = [5,0,10];
Z = V(2+sign(X))
Z = 1×9
0 5 0 10 10 5 10 5 0
  1 Comment
NN
NN on 5 May 2021
Thank you very much, i used this expression in an optimisation problem and i got this error
An error occurred while running the simulation and the simulation was terminated
Caused by:
  • sign

Sign in to comment.


Image Analyst
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
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.

Sign in to comment.


Walter Roberson
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
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;

Sign in to comment.


Matt J
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
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..

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by