I am currently working on an optimization problem to maximize the revenue from a system consisting of a PV plant. a battery, and an electrolyzer. When I try to solve it, I use "Intlinprog" but the code does not give any reslu*ts and the only thing to do is to force the code to stop ( Ctrl+C). I think the problems is in the constraints, becuase I used "Intlinprog" to solve a different problem and in that case it worked properly.
I also added to "Intlinprog" the option 'MaxTime' but in that case the results present a Relative Gap to high, I can not consder them as the optimal solutions.
Any assistance or insight on how to run "Intlinprog" properly would be greatly appreciated.
nome_simulazione='off_grid';
%% Import PRODUZIONE PV
filename = strcat(pwd , '\PRODUZIONE.csv');
delimiter = {''};
formatSpec = '%f%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'TextType', 'string', 'EmptyValue', NaN, 'ReturnOnError', false);
fclose(fileID);
PRODUZIONE = table(dataArray{1:end-1}, 'VariableNames', {'VarName1'});
clearvars filename delimiter formatSpec fileID dataArray ans;
Ppv = table2array(PRODUZIONE)';
% inserire degradazione annua 0.5%
Ppv1 = [Ppv Ppv*0.995 Ppv*0.99 Ppv*0.985 Ppv*0.980 Ppv*0.975 Ppv*0.97 Ppv*0.965 Ppv*0.960 Ppv*0.955 Ppv*0.95 Ppv*0.945 Ppv*0.94 Ppv*0.935 Ppv*0.93 Ppv*0.925 Ppv*0.920 Ppv*0.915]/1000000; % porto in [MW] da
%% Import PARAMETRI PER DEGRADAZIONE (la tecnologia usata è NMC)
load('aging_param_NMC.mat');
%% TEST VALUES:
PPVp = 100; % potenza nominale di picco del PV [MW]
HHV = 0.0394; % [MWh] potere calorifico di 1 kg di H2
dT = 0.25; % periodo di tempo tra un 'campione' ed il successivo [h]
incremento = 15; % incremento del costo per l'acquisto rispetto al prezzo di vendita
eta_pvg = 0.985*0.995; % rendimento scambio PV-Grid (inverter*traffo)
eta_pvb = 0.985*0.985; % rendimento scambio PV-Battery (inverter*inverter)
eta_be = 0.94; % rendimento scambio Battery-elettrolizzatore (inverter*traffo)
eta_pve=0.94; % rendimento di scambio PV-elettrolizzatore
Pemax=40;
Pbmax=40;
Cnom=8*Pbmax;
Cmin=75;
PriceH2=4;
startCostHIGH =800;
%startCostLOW = 0.2*Pemax;
SOCmin = 10; % valore di SOC minimo consentito
SOCmax = 90; % valore di SOC massimo consentito
SOCini = 50; % SOC iniziale
Price_electricity = 1;
% INITIALIZATIONS:
RicaviGiornalieri = 0;
ore_funz=0;
ore_funz_fine_vita=80000;
tot_day=0;
d = 0;
sum_fd = 0;
sum_fd_cy = 0;
sum_fd_cal = 0;
L = 0;
C_retention = Cnom;
Cnew = Cnom;
a= 0.008439;
b=-0.03224;
c=0.04967;
e=-0.07391;
f=0.997;
eta_battery = a*(Pbmax/Cnom)^4+b*(Pbmax/Cnom)^3+c*(Pbmax/Cnom)^2+e*(Pbmax/Cnom)+f;
eta_c = eta_battery; % rendimento di carica scelto in funzione del Crate massimo
eta_d = eta_c;
%tecnologia NMC (dati per lookup table modello simulink)
RTE = [0.00803226 0.011567 0.015622 0.025834]; % RTE
BP = [0.2 0.3 0.5 1]; % Breakpoints
%% OPTIMIZATION
while ore_funz <= ore_funz_fine_vita && (Cnew/Cnom)*100 >= Cmin
disp('day:')
disp(d+1)
daily_Ppv =(Ppv1(1,(d*96)+1:(d+1)*96)); % [MW] Potenza prodotta dal parco fotovoltaico nel giorno considerato
n = numel(daily_Ppv);
% rendimento elettrolizzatore
yrend = zeros(n,5);
yrend(:,1)=0;
yrend(:,2)=0.66;
yrend(:,3)=0.70;
yrend(:,4)=0.75;
yrend(:,5)=0.72;
% Potenza elettrolizzatore
yPelet=zeros(n,5);
yPelet(:,1)=0.04*Pemax;
yPelet(:,2)=0.25*Pemax;
yPelet(:,3)=0.5*Pemax;
yPelet(:,4)=0.75*Pemax;
yPelet(:,5)=1*Pemax;
% OPTIMIZATION PROBLEM
PowerFlow = optimproblem;
% VARIABLES
Ppvg = optimvar('Ppvg',1,n,'LowerBound',0,'UpperBound',daily_Ppv); % potenza PV -> rete
Ppve = optimvar('Ppve',1,n,'LowerBound',0,'UpperBound',Pemax/eta_pve); % potenza PV -> elettrolizzatore
Ppvb = optimvar('Ppvb',1,n,'LowerBound',0,'UpperBound',Pbmax/eta_pvb);
SOC = optimvar('SOC',1,n,'LowerBound',SOCmin,'UpperBound',SOCmax);
Pbe = optimvar('Pbe',1,n,'LowerBound',0,'UpperBound',Pbmax/eta_d);
Peletlogical = optimvar('Peletlogical',n,11,'Type','integer','LowerBound',0,'UpperBound',1);
z = optimvar('z',n,'Type','integer','LowerBound',0,'UpperBound',1);
Z_X = optimvar('Z_X',1,n,'Type','integer','LowerBound',0,'UpperBound',1); % charging binary variable array
Z_Y = optimvar('Z_Y',1,n,'Type','integer','LowerBound',0,'UpperBound',1); % discharging binary variable array
% OBJECTIVE
RPVG = sum(((Ppvg.*eta_pvg).*Price_electricity)*dT);
RH2 = sum((((Peletlogical(:,1).*yPelet(:,1)).*yrend(:,1))+...
((Peletlogical(:,2).*yPelet(:,2)).*yrend(:,2))+...
((Peletlogical(:,3).*yPelet(:,3)).*yrend(:,3))+...
((Peletlogical(:,4).*yPelet(:,4)).*yrend(:,4))+...
((Peletlogical(:,5).*yPelet(:,5)).*yrend(:,5)))/HHV)*dT*PriceH2;
startingCost = sum((z.*startCostHIGH)*dT);
PowerFlow.Objective = +RPVG-RH2+startingCost;
w = optimexpr(n,1);
idx =2:n;
if d==0
w(1,1)=Peletlogical(1,1);
w(idx,1)=Peletlogical(idx,1)-Peletlogical(idx-1,1)+...
Peletlogical(idx,2)-Peletlogical(idx-1,2)+...
Peletlogical(idx,3)-Peletlogical(idx-1,3)+...
Peletlogical(idx,4)-Peletlogical(idx-1,4)+...
Peletlogical(idx,5)-Peletlogical(idx-1,5);
end
if d>=1
w(1,1)= Peletlogical(1,1)-Results{1,d}.Peletlogical(n,1)+...
Peletlogical(1,2)-Results{1,d}.Peletlogical(n,2)+...
Peletlogical(1,3)-Results{1,d}.Peletlogical(n,3)+...
Peletlogical(1,4)-Results{1,d}.Peletlogical(n,4)+...
Peletlogical(1,5)-Results{1,d}.Peletlogical(n,5);
w(idx,1)=Peletlogical(idx,1)-Peletlogical(idx-1,1)+...
Peletlogical(idx,2)-Peletlogical(idx-1,2)+...
Peletlogical(idx,3)-Peletlogical(idx-1,3)+...
Peletlogical(idx,4)-Peletlogical(idx-1,4)+...
Peletlogical(idx,5)-Peletlogical(idx-1,5);
end
switchcons = w - z <= 0;
% CONSTRAINTS
PowerFlow.Constraints.Ch_on = Ppvb <= (Pbmax/eta_pvb)*Z_X;
PowerFlow.Constraints.Dis_on = Pbe <= (Pbmax/eta_d)*Z_Y;
PowerFlow.Constraints.notPositive = Z_X + Z_Y <= 1;
PowerFlow.Constraints.switchcons = switchcons;
PowerFlow.Constraints.tot=Peletlogical(:,1)+...
Peletlogical(:,2)+...
Peletlogical(:,3)+...
Peletlogical(:,4)+...
Peletlogical(:,5)<=1;
PowerFlow.Constraints.tot2=((Peletlogical(:,1).*yPelet(:,1)))+...
((Peletlogical(:,2).*yPelet(:,2)))+...
((Peletlogical(:,3).*yPelet(:,3)))+...
((Peletlogical(:,4).*yPelet(:,4)))+...
((Peletlogical(:,5).*yPelet(:,5)))-Pbe'== +Ppve';
PowerFlow.Constraints.PPV = daily_Ppv == Ppvb + Ppvg + Ppve;
PowerFlow.Constraints.rampa=z+...
Peletlogical(:,1)*2+...
Peletlogical(:,2)*2+...
Peletlogical(:,3)+...
Peletlogical(:,4)*2+...
Peletlogical(:,5)*2<=2;
PowerFlow.Constraints.SOC = optimconstr(n); % vincoli sullo Stato di Carica
for j=1:n
if j==1
PowerFlow.Constraints.SOC = SOC(j) == SOCini;
else
PowerFlow.Constraints.SOC(j) = SOC(j) == SOC(j-1)+((Ppvb(j)*eta_pvb)-(Pbe(j)*eta_d))/Cnew*(100*dT);
end
end
PowerFlow.Constraints.SOC2= SOC(end)==50;
options = optimoptions('intlinprog','Display','final','MaxTime',7);
[Result,fval,exitflag,output]= solve(PowerFlow,'options',options);
% RESULTS SAVING
SOCini = Result.SOC(length(SOC));
Results(d+1) = {Result}; % salvataggio dei profili di potenza giornalieri ottenuti dall'ott.
RicaviGiornalieri(d+1) = -fval;
index = zeros(1,length(Result.Peletlogical));
index(Result.Peletlogical >= 0.1)=1;
ore_funz = ore_funz+sum(index)/4;
tot_day=tot_day+length(daily_Ppv)/96;
years_of_life = (tot_day)/365;
d=d+1;
% BATTERY DEGRADATION:
% BATTERY TEMPERATURE EVALUATION (Simulink)
t = 0:0.25:23.75; % vettore dei tempi per il giorno
P_batt = timeseries((Result.Pbe-Result.Ppvb),t*3600); % timeseries per Simulink [MW]
T_amb = 22; % [°C]
%dati per tecnologia NMC
n_module = Cnom/0.013; % numero di moduli utilizzati
s = (1.162*0.110*2)+(0.445*0.110*2); % superficie di scambio per ogni modulo [m^2]
cell_area = n_module*s; % [m^2]
h_conv = 200; % heat transfer coefficient (convezione forzata)
w = 89; % [kg] peso cella
cell_mass = n_module*w;
cell_Cp_heat = 1600; % specific heat [J/kg/K]
tic
t_simulation = t(length(t))*3600;
if d==0
T_init = T_amb+273.15;
end
if d>=1
T_init = T_bat.Data(end)+273.15; % in ogni giorno T ini. = T fin. giorno precedente
end
sim('modello_Temperatura_Batteria')
toc
% fattori di stress:
% tecnologia NMC
[fd_aging_cy,fd_aging_cal,S_temp,Temp] = battery_degradation(Result,aging_param_NMC,T_bat,t);
sum_fd_cy = sum_fd_cy+fd_aging_cy;
sum_fd_cal = sum_fd_cal+fd_aging_cal;
sum_fd = sum_fd+fd_aging_cy+fd_aging_cal;
L(d+1) = 1-aging_param_NMC.alfasei*exp(-aging_param_NMC.betasei*sum_fd)-(1-aging_param_NMC.alfasei)*exp(-sum_fd);
Cnew = Cnom*(1-L(d+1)); % nuovo valore di capacità per il giorno successivo
C_retention(1,d+2) = Cnew; % profilo di variazione della capacità
disp('nuovo valore di capacità:')
disp(Cnew)
end
save('off_grid_C/D')
