Energy storage optimisation problem - separate charge and discharge constraint

Hi all!

I am currently working on an optimization problem to maximize the revenue from a combined wind turbine and energy storage system. With the code below, the system charges and discharges simultaneously at certain times. I am able to change the charge and discharge variables to a single variable with the lowerbound as a negative and positive number corresponding to discharging and charging, respectively. But in doing this, I cannot accurately apply the charge and discharge efficiencies of the storage device.

prob = optimproblem;

% Decision variables

% Energy storage system decision variables
ESS_ch = optimvar('ESS_ch',T,'LowerBound',0,'UpperBound',ESS_Pmax);
ESS_disch = optimvar('ESS_disch',T,'LowerBound',0,'UpperBound',ESS_Pmax);
ESS_SOC = optimvar('ESS_SOC',T,'LowerBound',0,'UpperBound',ESS_Cmax);

% Output power to the grid variables
Grid_E = optimvar('Grid_E',T,'LowerBound',0);

% Energy storage operational constraints
prob.Constraints.energyStorage  = optimconstr(T);
prob.Constraints.energyStorage(1) = ESS_SOC(1) == 10;
prob.Constraints.energyStorage(2:T) = ESS_SOC(2:T) == ESS_SOC(1:T-1)*(1-d) - ESS_disch(1:T-1)/eff + ESS_ch(1:T-1)*eff;

% Energy balance constraint (WT_E = wind turbine power output as vector of
% length T)
prob.Constraints.EnergyBalance = Grid_E == WT_E + ESS_disch - ESS_ch;

prob.ObjectiveSense = 'maximize';
prob.Objective = sum(Grid_E.*Price);

I have also attempted to apply a constraint that dictates the charge power multilied by the discharge power of the storage device is equal to zero at any single time step, see below.

prob.Constraints.chargeonoff  = optimconstr(T);
for i = 1:T
    prob.Constraints.chargeonoff(i) = ESS_disch(i)*ESS_ch(i) == 0;
end

However, when I add this the following error occurs.

Error using optim.problemdef.OptimizationProblem/solve

SOLVE requires a non-empty initial point structure to solve a nonlinear problem.

Any assistance or insight on how to force the program to either charge or discharge at each time step would be greatly appreciated.

[Shehbaz Aslam]

Hope you are fine I need help on optimization toolbox. [e-mail]

Thanks

[James Cresswell]

Hi Shehbaz,

I am not especially experienced with the optimization toolbox and still figuring things out for myself, so unfortunately I don't think I will be able to help.

Best,

James

[NN]

Hi James,

I am working on microgrid optimization with solar and energy storage.I am also checking the Matlab example file Microgrid EMS optimization as a reference

https://www.mathworks.com/matlabcentral/fileexchange/73139-microgrid-energy-management-system-ems-using-optimization?s_tid=srchtitle

There are some queries which I didn't understand. I request you all to help me to understand the code and then only I will be able to do my model. These are a few queries which I like to ask.

In the Matlab example, they load a file PvLoadPrice.I understand the load details are for variable load given in the system.PV has clear and cloudy data. Can I know how do we use such data?in which format? From which website we download it for research purposes?Also, how do we separate such data against each parameter? How do I forecast those data, say for the price how do I compare price details before and after a load variation(i have load fluctuation event in my model ) ? should I do forecasting of the price? I want to minimize the price of grid net exchange.plz advice. Can I have any resources to understand the code written in this example file?

The energy storage charges and discharges as per the power availability in the grid. So how can I model the optimization problem? and what should be the constraints? I have already modeled a microgrid with charging and discharging energy storage mechanism as per the load availability. So how can I link the model with the optimization code?

I request James, ALi and Alan Weiss to look into thsi matter and help .

Thanks in advance.

[Joan Ernesto]

How did you manage to make it continuous? I'm trying to do something similar to minimize the grid cost, but at some timesteps the optimizer decides to discharge the battery at random timestep instead of continuously discharging it.

Answer 1

[Alan Weiss]

The error message is spot on: when you multiply variables, the problem becomes nonlinear. However, I am sure that you would rather keep the problem linear if you can.

Choose a value M that upper bounds a nonnegative decision variable x. I mean x is ESS_disch or ESS_ch. Create a binary variable z_x that will indicate when x is positive. Set a linear constraint so that z tracks whether x is positive or not.

x <= M*z_x; % if x is positive, z_x has to be positive

Include z in the objective function, so the objective is minimized when z_x = 0.

Now to make the constraint that x and y cannot both be positive, create the linear constraint

z_x + z_y <= 1;

I may have some detail wrong, but I think that you see the idea. You are creating extra integer variables and linear constraints, and turning the problem from a linear programming problem to a mixed-integer linear programming problem. That keeps the problem out of the general nonlinear framework, which is a good thing.

Alan Weiss

MATLAB mathematical toolbox documentation

[James Cresswell]

Hi Alan,

Thanks very much for your reply. Yes, I very much want to keep the program linear. I have attempted to apply your suggestions but I am getting a bit confused with some parts.

I created the decision variables for ESS_ch and ESS_disch.

M=30;
ESS_ch = optimvar('ESS_ch',T,'LowerBound',0,'UpperBound',M);
ESS_disch = optimvar('ESS_disch',T,'LowerBound',0,'UpperBound',M);

Next I created two binary variable (one for ESS_ch and one for ESS_disch) which I assume was right, or is only one binary variable necesary? And is it meant to be a 1x1 or Tx1 variable? So I'm a little unsure about the constraints as well.

%create binary variables   *************************
Z_X = optimvar('Z_X',T,'Type','integer','LowerBound',0,'UpperBound',1); % charging binary variable array  
Z_Y = optimvar('Z_Y',T,'Type','integer','LowerBound',0,'UpperBound',1); % discharging binary variable array 

% Consraints for binary variables
prob.Constraints.Ch_on = ESS_ch <= M*Z_X;
prob.Constraints.Dis_on = ESS_disch <= M*Z_Y;
prob.Constraints.notPositive = Z_X + Z_Y <= 1;
% ***************************************************

Finally, I do not know where to insert the binary variables in order to maximize the objective, I get different error messages depending on where I try to put them. Could you explain a bit more please? I thought the binary variables would have needed to be included in the energyStorage SOC constraint as well as the objective, but you also said that I shouldn't multiply a variable by a variable.

prob.Constraints.energyStorage(2:T) = ESS_SOC(2:T) == ESS_SOC(1:T-1)*(1-d) + ESS_ch(1:T-1) - ESS_disch(1:T-1);

prob.Constraints.EnergyBalance = Grid_E == WT_E - Z_X.*ESS_ch + Z_Y.*ESS_disch;

prob.ObjectiveSense = 'maximize';
prob.Objective = sum(Grid_E.*Price); % Price = Tx1 vector

Many thanks in advance

[Alan Weiss]

You got it right, Z_X and Z_Y are vectors the same sizes as Ess_ch and Ess_disch respectively. You don't need to include these variables in any constraints except in the first three constraints:

prob.Constraints.Ch_on = ESS_ch <= M*Z_X;
prob.Constraints.Dis_on = ESS_disch <= M*Z_Y;
prob.Constraints.notPositive = Z_X + Z_Y <= 1;

They are not needed in energyStorage or EnergyBalance, I believe. You could put them in the objectiive like this if you want:

prob.Objective = sum(Grid_E.*Price) + sum(Z_X) + sum(Z_Y); % Price = Tx1 vector

Again, I am not 100% sure that this formulation is correct, but I think that it is worth a try.

Alan Weiss

MATLAB mathematical toolbox documentation

[James Cresswell]

Great, that works now!

Thanks very much for your help. :)

James

[Alan Weiss]

I'm happy that my suggestion worked, and that you let us know. I was a little uncertain that I had understood everything correctly, so it is especially gratifying to hear that it all worked for you.

I think that I did say something wrong. Because you are maximizing, the objective function should probably be

prob.Objective = sum(Grid_E.*Price) - sum(Z_X) - sum(Z_Y); % Price = Tx1 vector

because that tries to make Z_X and Z_Y equal to zero, when possible. But it probably doesn't matter for your problem.

Alan Weiss

MATLAB mathematical toolbox documentation

[Ali Almarzooqi]

Dear James,

I am working on a similar optimization problem except that my renewable source is Solar. I wrote my code in a different way, but it’s still not working correctly.

I would like to use your code if you don’t mind, and can I know what “d” represent in the SOC constraint?

Best regards,

Ali

[James Cresswell]

Hi Ali,

I do not mind, I got a lot of inspiration from the matlab community example below anyway. It uses solar PV input so it might be quite applicable for you. In the example however, it assumes an ideal systme by representing charge/discharge as a single variable with positive and negative upper and lower bounds, respectively. Which doesn't allow for applying the charge/discharge efficiency.

https://se.mathworks.com/matlabcentral/fileexchange/73139-microgrid-energy-management-system-ems-using-optimization?s_tid=srchtitle

"d" is a constant in my code for the standing losses (self-discharge) of the energy storage. Typically, this is <=5% per month for Li batteries, and then scale that down to what that is per hour or half-hour period of your system.

Best,

James

[Ali Almarzooqi]

Hi James,

Thank you a lot for your clarification and support!! Appreciated.

Best regards,

Ali

[NN]

Dear Sir,

Below is the link for an example file from matlab,

https://www.mathworks.com/matlabcentral/fileexchange/73139-microgrid-energy-management-system-ems-using-optimization?s_tid=srchtitle.

As if in the link, I am trying to do energy management system for EV batteries of 100 vehicles.I am considering each vehicle 100kWh and 20kW.In the example for single battery they used the below code for energy contraint,

% Battery SOC Energy constraints (keep between 20%-80% SOC)

battEnergy = 3.6e6*BattCap;

batteryMinMax.Emax = 0.8*battEnergy;

batteryMinMax.Emin = 0.2*battEnergy;

for my case considering discharge and charging batteries can I use the below code:

EVregBattCap = 100;

EVchBattCap = 100;

EVregbatteryMinMax.Pmin = -20e3;

EVregbatteryMinMax.Pmax = 20e3;

EVdschbattEnergy = 3.6e6*EVdchBattCap;

EVdschbatteryMinMax.Emax = 0.8*EVdschbattEnergy;

EVdschbatteryMinMax.Emin = 0.2*EVdschbattEnergy;

EVchbattEnergy = 3.6e6*EVchBattCap;

EVchbatteryMinMax.Emax = 0.8*EVchbattEnergy;

EVchbatteryMinMax.Emin = 0.2*EVchbattEnergy;

Or should I multiply it with 100 no:s?

When I run the system, I am expected to receive optimization signal as +20 e3 /0/-20e3. but I get a higher value.

So as in the example file, instead of a single battery when I consider multiple batteries, please advise the changes I should incorporate.

[Alan Weiss]

I am not sure, but you might need to set constraints on entire vectors, not on scalar values. For example, I suppose that you have a vector of battery energy, and you need to set the constraints (bounds, rate of flow) for each battery, meaning a vector of constraints. But without seeing your code I cannot be sure how you impose the constraints, so I might be wrong.

Alan Weiss

MATLAB mathematical toolbox documentation