22 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

Matt J
on 5 May 2021

Edited: Matt J
on 5 May 2021

Replace

sqrt(x.^2)./x

with

sign(x)

Walter Roberson
on 5 May 2021

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

Walter Roberson
on 6 May 2021

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?

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

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

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

Start Hunting!