How can I implement an "if" in optimization problem?

In the example https://www.mathworks.com/matlabcentral/fileexchange/73139-microgrid-energy-management-system-ems-using-optimization?s_tid=srchtitle energy management optimization has been done in presence of an energy storage system. The optimproblem has been set as follow:

%% Define 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);

%% Objective Function

% Minimize cost of electricity from the grid

prob.ObjectiveSense = 'minimize';

% prob.Objective = dt*C'*PgridV - FinalWeight*EbattV(N);

prob.Objective = dt*C'*PgridV - FinalWeight*EbattV(N);

%% Constraints

% 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;

Last line describes energy storage behaviour, without considering the efficiency. Setting eta_c and eta_d as charge and discharge efficiency rispectively, how can I implement an "if" problem in such a constraint? I mean something like:

if PbattV>0 prob.Constraints.energyBalance(2:N) = EbattV(2:N) == EbattV(1:N-1) - PbattV(1:N-1)*dt/eta_d;

if PbattV<0 prob.Constraints.energyBalance(2:N) = EbattV(2:N) == EbattV(1:N-1) - PbattV(1:N-1)*dt*eta_c;

I don't manage to figure it out because PbattV is an optimvar.

Please, do you have any advice?

Answer 1

[Alan Weiss]

Yoou might be able to find some techniques in the new documentation topic Integer and Logical Modeling.

Alan Weiss

MATLAB mathematical toolbox documentation

[Danny]

Hi Alan,

thanks for your advice and attention. I've followed your suggestion and I have considered the battery as reservoir flow problem formulated in the documentation. I've set the problem as follow:

% Decision variables

PgridV=optimvar('PgridV',N);

EbattV = optimvar('EbattV',N,'LowerBound',batteryMinMax.Emin,'UpperBound',batteryMinMax.Emax);

M=400e3; %max power rate

Pbatt_ch = optimvar('Pbatt_ch',N,'LowerBound',0,'UpperBound',M);

Pbatt_disch = optimvar('Pbatt_disch',N,'LowerBound',0,'UpperBound',M);

%Binary variables

z_in = optimvar('z_in',N,'Type','integer','LowerBound',0,'UpperBound',1); % charging binary variable array

z_out= optimvar('z_out',N,'Type','integer','LowerBound',0,'UpperBound',1); % discharging binary variable array

%% Constraints

%Binary variables

prob.Constraints.charge_on = Pbatt_ch <= M*zin_;

prob.Constraints.discharge_on = Pbatt_disch <= M*z_out;

prob.Constraints.notbothpositive = z_in + z_out <= 1;

% 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) +(Pbatt_ch(1:N-1)*eta_c)*dt-(Pbatt_disch(1:N-1)/eta_d)*dt;

%%Objective function (minimize cost of electricity from the grid)

prob.ObjectiveSense='minimize';

prob.Objective=dt*C'*PgridV-FinalWeight*EbattV(N);

I'm not pretty sure about last line of energyBalance constraint and I don't know if I got it.

Best Regards,

Danny