I am getting non-linear constraint error in the optimization problem

I am getting this error while solving the optimization problem using 'intlinprog' solver. The error I am getting is:

Error using optim.internal.problemdef.ProblemImpl/prob2structImpl

Problems with integer variables and nonlinear equality

constraints are not supported.

Error in optim.internal.problemdef.ProblemImpl/solveImpl

Error in optim.problemdef.OptimizationProblem/solve

Error in Test3_some_equality_removed (line 193)

[min_lcohsol,fval,exitflag,output] = solve(min_lcoh,'options',options);

The code I am trying to run is:

% number of time steps
nPeriods = length(GHI);

%% PV, Electrolyser and Battery Min Max values
% Pv
Ppv_rated_min = 0;
Ppv_rated_max = 800; % MW

% Electrolyser
Pel_rated_min = 0;
Pel_rated_max = 500; % MW

% Battery
ES_rated_min = 0;
ES_rated_max = 500; % MWh

%% Introducing decision variables
del = optimvar('del',nPeriods,1,'Type','integer','LowerBound',0,'UpperBound',1);
P_pv = optimvar('P_pv',nPeriods,1,'Type','continuous','LowerBound',0,'UpperBound',Ppv_rated_max);
P_el = optimvar('P_el',nPeriods,1,'Type','continuous','LowerBound',0,'UpperBound',Pel_rated_max);
P_ch = optimvar('P_ch',nPeriods,1,'Type','continuous','LowerBound',0,'UpperBound',ES_rated_max);
P_dch = optimvar('P_dch',nPeriods,1,'Type','continuous','LowerBound',0,'UpperBound',ES_rated_max);
Ppv_rated = optimvar('Ppv_rated',nPeriods,1,'Type','integer','LowerBound',Ppv_rated_min,'UpperBound',Ppv_rated_max);
Pel_rated = optimvar('Pel_rated',nPeriods,1,'Type','integer','LowerBound',Pel_rated_min,'UpperBound',Pel_rated_max);
ES_rated = optimvar('ES_rated',nPeriods,1,'Type','integer','LowerBound',ES_rated_min,'UpperBound',ES_rated_max);


%% Power Balance
Pb_ineq = P_pv >= P_el + P_ch - P_dch;

%% PV, Electrolyser and battery sizing inequality
% PV
PV_rating1 = Ppv_rated_min <= Ppv_rated;
PV_rating2 = Ppv_rated <= Ppv_rated_max;

% Electrolyser
Electrolyser_rating1 = Pel_rated_min <= Pel_rated;
Electrolyser_rating2 = Pel_rated <= Pel_rated_max;

% Battery inequality
Battery_rating1 = ES_rated_min <= ES_rated;
Battery_rating2 = ES_rated <= ES_rated_max;


%% Constraints particular to PV
P_pv = f_pv.*(Ppv_rated.*1e6).*(GHI./(G_std.*1000)).*(1 + alpha_pv.*(T_cell-T_cell_std))/(1e6); % MW


%% Constraints particular to Electrolyser
% initializing size of mass of hydrogen vector
mh2 = zeros(nPeriods,1);

Pel_rated_aux = zeros(nPeriods,1);

% introducting binary variable to define on/off state
Pel_eq1 = Pel_rated_aux == Pel_rated.*del;

% inequalities required after introducing binary variable
el1 = Pel_rated_aux <= Pel_rated - (1-del).*Pel_rated_min;
el2 = Pel_rated_aux >= Pel_rated - (1-del).*Pel_rated_max;
el3 = Pel_rated_aux <= Pel_rated_max.*del;
el4 = Pel_rated_aux >= Pel_rated_min.*del;

% for electrolyser modulation
Electrolyser_modulation1 = (Pel_min.*Pel_rated_aux) <= P_el;
Electrolyser_modulation2 = P_el <= (Pel_max.*Pel_rated_aux);

%% Constriants Particular to Hydrogen Production
mh2 = (nel.*P_el.*1000)./lhv;
mh2_eq2 = sum(mh2) >= 400;


%% Constraints particular to Battery
% Initial Charge
SOC_in = 0.5;
ES = zeros(nPeriods,1);

% for Battery Modulation
Battery_modulation1 = (ESmin.*ES_rated)<= ES;
Battery_modulation2 = ES <= (ESmax.*ES_rated);

% State of the charge in battery
ESeq1 = ES(1) == ES_rated.*SOC_in;
ESeq2 = ES(2:nPeriods) == ES(1:nPeriods-1).*(1-sigmaES) + nESch.*P_ch(2:nPeriods) - P_dch(2:nPeriods)./nESds;

%% Economic Part
%Electrolyser
year = 2024;
ko = 301.04; k = 11603; a = 0.649; b = -27.33; 
vo = 2020; v=year;

% electrolyzer capex constraints
q = Pel_rated(1).*1000; % in kW
ELectrolyser_capex = 1.*(ko + (k./q).*(q.^a)).*(v./vo).^b; 
Electrolyser_opex = 0.075.*ELectrolyser_capex;

% PV 
PV_capex = pv_per_kw.*Ppv_rated(1);
PV_opex = 0.01.* PV_capex;

% Battery
Battery_capex = battery_per_kwh.*ES_rated(1);
Battery_opex = 0.025.*Battery_capex;

% Burner Replacement
Burner_power = (nel.*Pel_rated(1)); % Peak hydrogen production per hour multiplied by lhv,MW
Burner_capex = (0.1896 * Burner_power.^0.6113 - 0.001308).*1.09; % included euro to dollor conversion

% System installation
System_installation_capex = 0.12.*(ELectrolyser_capex+PV_capex+Battery_capex+Burner_capex);

%Overall Capex
Capex = ELectrolyser_capex+PV_capex+Battery_capex+Burner_capex+System_installation_capex;

%Overall Opex
Opex = Electrolyser_opex+PV_opex+Battery_opex;

% Lcoh
% Computing the two sums without loops
one = sum(Opex ./ discount_factor);
two = sum((sum(mh2) .* ((1 - Stack_dr).^n)) ./ discount_factor);

% Computing LCOH
lcoh = (Capex + one) / two;

%% Optimization Process
min_lcoh = optimproblem('ObjectiveSense','minimize');

% Objective
min_lcoh.Objective = lcoh;

% Constraints
min_lcoh.Constraints.Pb_ineq = Pb_ineq;
min_lcoh.Constraints.PV_rating1 = PV_rating1;
min_lcoh.Constraints.PV_rating2 = PV_rating2;
min_lcoh.Constraints.Electrolyser_rating1 = Electrolyser_rating1;
min_lcoh.Constraints.Electrolyser_rating2 = Electrolyser_rating2;
min_lcoh.Constraints.Battery_rating1 = Battery_rating1;
min_lcoh.Constraints.Battery_rating2 = Battery_rating2;
min_lcoh.Constraints.Pel_eq1 = Pel_eq1;
min_lcoh.Constraints.el1 = el1;
min_lcoh.Constraints.el2 = el2;
min_lcoh.Constraints.el3 = el3;
min_lcoh.Constraints.el4 = el4;
min_lcoh.Constraints.Electrolyser_modulation1 = Electrolyser_modulation1;
min_lcoh.Constraints.Electrolyser_modulation2 = Electrolyser_modulation2;
min_lcoh.Constraints.mh2_eq2 = mh2_eq2;
min_lcoh.Constraints.Battery_modulation1 = Battery_modulation1;
min_lcoh.Constraints.Battery_modulation2 = Battery_modulation2;
min_lcoh.Constraints.ESeq1 = ESeq1;
min_lcoh.Constraints.ESeq2 = ESeq2;

options = optimoptions('intlinprog','Display','final');
[min_lcohsol,fval,exitflag,output] = solve(min_lcoh,'options',options);

Please help me solve this problem. I am not able to figure out where the problem actually is.

[Sam Chak]

Hi @Anup, I was wondering whether the displayed error messages are sufficiently informative for both yourself and other ordinary MATLAB users when troubleshooting and fixing the code. While your code appears to be relatively organized, it would greatly benefit from detailed annotations describing the purpose of each line.

Answer 1

[Torsten]

"solve" chooses MATLAB's "ga" to tackle your problem, but "ga" has the restriction that integer constraints and nonlinear inequaltity constraints cannot be imposed simultanously.

From the documentation of "ga":

Note

When there are integer constraints, ga does not accept nonlinear equality constraints, only nonlinear inequality constraints.

Thus you will somehow have to reformulate your problem in order to solve it using MATLAB.