2 views (last 30 days)

Show older comments

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

Alan Weiss
on 2 Jun 2021

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

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

Start Hunting!