Problem-Based Optimization, No feasible solution found

In this simplified version of a system designed to produce hydrogen from pv, I want to find out the lowest possible capacity of the electrolyzer?

The problem is that when I run the model it says that there is no solution avaliable. However, when I fill in the capacity of the electrolyzer by hand I can easily find a solution for the constraints.

%% Problem-Based Optimization 
load('CF_pv.mat');                  % Hourly data of capacity factor sun

energyprob = optimproblem;

PV = 748;                                                                               % in MW
Electrolyzer = optimvar('Electrolyzer',1,'LowerBound',0);                               % Capacity of ELectrolyzer in MW

Electrolyzer_costs =[(200+0.7*200)];
costs = Electrolyzer_costs*Electrolyzer;

%Objective Function
energyprob.Objective = costs;

%% Optimization constraints
PV_Load = PV*CF_pv;

T = 8760;                                                      % hours in a year
Bat_In = optimvar('Bat_In',T,'LowerBound',0);                  % Electricity send to battery     
Bat_Out = optimvar('Bat_Out',T,'LowerBound',0);                % Electricity send from battery to Electrolyzer 

Bat_In_c = optimconstr(T);
Bat_Out_c = optimconstr(T);

for t=1:T
    Bat_In_c(t) = Bat_In(t) == (PV_Load(t)-1.1*Electrolyzer);       % Abundant Electricity, above 110% of the limit of the electrolyzer
    Bat_Out_c(t) = Bat_Out(t) == 0.1*Electrolyzer-PV_Load(t);       % Electricity necessary to keep electrolyzer running    
end

Elec_In = PV_Load - Bat_In + Bat_Out;                               % Generated electricity - input in battery + electricty from battery to electrolyzer     

energyprob.Constraints.Bat_In_c = Bat_In_c;
energyprob.Constraints.Bat_Out_c = Bat_Out_c;

Desired_Elec_In = 1.7618*10^6;                                      % Desired electricty input for a whole year

cons1 = Elec_In == sum(Desired_Elec_In);

energyprob.Constraints.cons1 = cons1;

[sol,fval] = solve(energyprob);

sol.Electrolyzer % display electrolyzer capacity

Answer 1

[Alan Weiss]

Your problem is infeasible as given. Using the routines in Investigate Linear Infeasibilities, I found the following. Convert the problem to LP matrices:

problem = prob2struct(energyprob);

Add a null Aineq and bineq inequality constraint.

Aineq = zeros(1,size(problem.Aeq,2));
bineq = 0; % The code needs an Aineq argument of the correct size

I then ran the code in Remove iIS In A Loop, and found that the following are conflicting lines in the problem.Aeq matrix together with the problem.beq data (conflict with the bounds):

Row 8760 and 17513

Row 8759 and 17512

Row 8758 and 17511

etc.

It just kept spitting them out, so I stopped the process. Armed with this knowledge I am sure that you can find the error in your setup.

By the way, I recently found some errors in the published code for the generatecover function. The corrections you must make in the generatecover code:

% Look for this line
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq),beq(activeAeq),lb,ub);
% change it to
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);

% Look for this line
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
% change it to
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);

But I suggest that you do not run generatecover on this problem; I have been running it for hours and it is still not finished.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Alan Weiss]

Just a bit more info. Looking at Rows 8760 and 17513, each has only two nonzero entries. 8760:

Aeq(8760,8760) = 1.0000

Aeq(8760,17521) = 1.1000

beq(8760) = 0

Together, these mean that x(8760) + 1.1*x(17521) = 0. But the lower bound on all x components is 0, so this implies x(8760) = x(17521) = 0.

Examine row 17513.

Aeq(17513,17513) = 1

Aeq(17513,17521) = -0.1

beq(17513) = -115.8937

These mean x(17513) - 0.1*x(17521) = -115.8937.

Therefore we must have x(17521) > 0 (actually x(17521) > 1158.937). But we already found that x(17521) = 0. So this is a contradiction.

Alan Weiss

MATLAB mathematical toolbox documentation

[Rogier Doodeman]

Thank you for your answer, this was indeed the problem.

For those interested the code delvired the right answer when I removed the contraints determining the Bat_In and Bat_Out value and instead wrote contraints 2 and 3.

%% Problem-Based Optimization 
load('CF_pv.mat');                  % Hourly data of capacity factor sun

energyprob = optimproblem('ObjectiveSense','minimize');

PV = optimvar('PV','Type','integer','LowerBound',500);                     % Capacity of pv panels in MW
Electrolyzer = optimvar('Electrolyzer','Type','integer','LowerBound',0);                               % Capacity of ELectrolyzer in MW

PV_costs =[(550+5.5)*10^3]*PV;                                  %[(550+5.5)*(max(Output_Solar_DC)] times the peak output so the price is per meter
Electrolyzer_costs =(200+0.7*200)*Electrolyzer;

%% Optimization constraints
PV_Load = PV*CF_pv;
Desired_Elec_In = 1.7618*10^6;                                      % Desired electricty input for a whole year

T = 8760;
M = 500000;         % Value that can not be surpassed
Bat_In = optimvar('Bat_In',T,'LowerBound',0,'UpperBound',M);
Bat_Out = optimvar('Bat_Out',T,'LowerBound',-M,'UpperBound',0);
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 

Bat_Mut = Bat_In + Bat_Out;                                 % Battery mutation per hour, Bat_Out are all negative values so it will subtract 
Elec_In = PV_Load - Bat_Mut;

%Objective Function
energyprob.Objective = PV_costs + Electrolyzer_costs +sum(Z_X) + sum(Z_Y); % + Battery_Ene_costs + Battery_Pow_costs

cons1 = sum(Elec_In) == Desired_Elec_In;
cons2 = Elec_In <= 1.1*Electrolyzer;
cons3 = Elec_In >= 0.1*Electrolyzer;

cons4 = Bat_In <= M*Z_X;                                    
cons5 = Bat_Out >= -M*Z_Y;                                  
cons6 = Z_X + Z_Y <= 1;                                     

energyprob.Constraints.cons1 = cons1;
energyprob.Constraints.cons2 = cons2;
energyprob.Constraints.cons3 = cons3;

energyprob.Constraints.cons4 = cons4;
energyprob.Constraints.cons5 = cons5;
energyprob.Constraints.cons6 = cons6;

[sol,fval] = solve(energyprob)
sol.PV 
sol.Electrolyzer % display electrolyzer capacity

To make sure the battery does not charge and discharge at the same time constraints 4 - 6 are used. For more information on this you can look at this thread. https://nl.mathworks.com/matlabcentral/answers/577642-energy-storage-optimisation-problem-separate-charge-and-discharge-constraint