how can i correct intlingprog error, when identifying an integer problem as non-linear?

Hi Matlab community,

i ran the following code in matlab, but the intlingprog identify the problem as non-linear, please how can i rectify it?

Thanks

the code:

%============= program for multi-energy system MES sizing=============
% this program consider optimal selection from the sub-categories of MES
% equipment, i.e CHP, PV, HP, GB. P2G, FUEL CELL
clc
% Load all the datas
LD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...
    4 2.5 2.3 0.1]; % Home Load Data for 24 hrs in kW
HD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...
    4 2.5 2.3 0.1]; % Heat Load Data for 24 hrs in kW
CD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...
    4 2.5 2.3 0.1]; % cooling---------- Load Data for 24 hrs in kW
SolarIrradiation=[0 0 0 0 0.1 0.3 0.5 0.65 0.8 0.9 1 0.95 0.9 0.8 0.65...
    0.4 0.3 0.08 0 0 0 0 0 0]; % for 24 hrs in kW/m^2
WindVelocity=repmat(5,1,24); % for 24 hrs in m/sec
Cgrid=[0.033 0.027 0.020 0.017 0.017 0.029 0.033 0.054 0.215...
    0.572 0.572 0.572 0.215 0.572 0.286 0.279 0.086 0.059 0.050 0.061...
    0.181 0.077 0.043 0.037]; % all prices and costs are in Euros
SellingCost=1.1; %feed in tarrif to the grid
Cgas=0.5;
% ==========PreDefined variables=================
%==========economy parameters
a=0.06;  %interest rate, residual value and proposed lifecycle
v=input('please input the carbon emission target'); % v is betwwen 0 to 1
Ta=input('please input the projected service life of MES equipment'); % Ta is 25years
Tloss=0.1; grideff=0.80;
AF=(a*((1+a)^Ta)/((1+a)^Ta)-1); %annuity factor

%====== parameters variables for ICE and Reciprocating combustion engine
gLHV=25.3; %low heating value of natural gas in kw/m3
RCEele=0.75; iceelec=0.65; %electrical efficiency of the two CHP
RCEloss=0.03; iceloss=0.04; %heat loss factor of the two CHP
rceyr=15; iceyr=20; %life span of each CHP
rcecost=200; icecost=250; %unit cost of investment of each CHP
rcemc=0.3; icemc=0.4; % proportion of unit cost for maintenance
Prcemin=5; Prcemax=15; Picemin=5; Picemax=15;
rcerep=(Ta/rceyr)-1; icerep=(Ta/iceyr)-1;

%======= parameters variables for PEMFC, SOFC and MOFC===============
hLHV=33.3; %low heating value of hydrogen gas in kw/m3
Pemeff= 0.7; sofeff=0.8; Mofceff= 0.88; %electrical efficiency of FC
phloss=0.3; sofhloss=0.4; mofhloss=0.5; %heat loss of considered FC
Pemmin=5; Pemmax=15; Psomin=5; Psomax=15; Pmomin=5; Pmomax=15;
pemyr=15; sofyr=20; mofyr=25; %life span of each FC
pemcost=250; sofcost=300; mofcost=350; %unit cost of investment
pemmc=0.03; sofmc=0.04; mofmc=0.05; % proportion of unit cost for maintenance
pemrep=(Ta/pemyr)-1; sofrep=(Ta/sofyr)-1; mofrep=(Ta/mofyr)-1;

% ==============parameters variables for electrolyser============
elez=10; %specific electrical energy demand of electrolyser
eleyr=15; %life span of electrolyser
elzcost=300; %unit cost of elz investment
elzmc=0.03; %proportion of unit cost for maintenance
Hemax=15; Hemin=5;
elzrep=(Ta/eleyr)-1;

% =============photovoltaic parameters=============
PR=270; % rated power of selected pv module
fR=0.8; %derating factor of pv module
TSC=25; %standard temperature of pv cell
GS=1000; %rated solar radiation at kw/m2
pvyr=25; %pv life span of pv
pvcost=200; %pv unit invesmtent cost
pvmc=0.05; %proportion of investment cost for annual maintenance
PVmin=5; PVmax=10; pvrep=(Ta/pvyr)-1;

% ===========Solar Thermal collector parameters============
TSeff=0.65;  %collector efficiency and installed area available
TSyr=25; TScost=200; TSmc=0.03; %collector cost, mcost and lifespan
TSAmin=0; TSAmax=300; TSrep=(Ta/TSyr)-1;

%============wind turbine parameters=============
Vci=4; % Cut-in Speed for wind turbine in m/sec
Vco=24; % Cut-out Speed for wind turbine in m/sec
Vr=14; % rated speed for wind turbine in m/sec
pw=5; %rated wind turbine power output
pwcost=250; %unit invetment cost of pv
pwmc=0.03; wtyr=20; %wind turbine mc cost and life span
wtrep=(Ta/wtyr)-1;

%===============parameters of heating equipment considered=================
Ashpeff=1.3; gshpef=1.4; %efficiency of air source & ground source hp
gbeff=0.92; gbheff=0.99; %efficiency of gas boiler and gas source hp
ashpyr=15; wshpyr=20; gbyr=15; gshpyr=15; %lifespan of each heat pump
ashcost=100; wshpcost=150; gbcost=200; gshpcost=250; %unit cost of hp
ashmc=0.03; wshpmc=0.04; gbmc=0.03; gshpmc=0.04; %maintenance cost
ashprep=(Ta/ashpyr)-1; gbrep=(Ta/gbyr)-1; wshprep=(Ta/wshpyr)-1;
gshphrep=(Ta/gshpyr)-1;

%===============parameters of cooling equipment considered=================
copAC=1.3; copEC=1.6; %coff. of performance of absorption and elec. chiller
ACyr=15; ECyr=20; %lifespan of chillers
ACcost=200; ECcost=250; %unit investment cost
ACmc=0.04; ECmc=0.03; %proportion of maintenance cost from investment cost 
ACrep=(Ta/ACyr)-1; ECrep=(Ta/ECyr)-1;

            % Energy storage variables
%======== electrical storage using electrochemical (battery)==============
echb=0.96; edchb=0.96;       %charging abd discharching efficiency of Battery
PbChmax=1;                  % Battery max allowed charging power in kW
PbDischmax=1;                % Battery max allowed discharging power in kW
selfdchb=0.001;                 %self-discharging parameters
BEcost=200; BEmc=0.03; BEyr=10; %BES unit cost, maintenance and life span
BEdgcost=0.1;                   %BES charging and discharging degradation cost
Brep=(Ta/BEyr)-1;

% ==========Thermal storage variables===========
qchmax=1;                   % thermal max allowed charging power in kW
qdchmax=1;                   % thermal max allowed discharging power in kW
qselfdch=0.005;             %self-discharging parameters
echq=1.00; edchq=1.00;      %charging abd discharching efficiency of Battery
tmin=10; tmax=35;            %thermal storage temperature
TEScost=200; TESmc=0.3; TESyr=10; %TES unit cost, maintenance and life span
TESdgcost=0.1;               %TES charging and discharging degradation cost
TESrep=(Ta/TESyr)-1;
  
% =============cold storage variables=============
SOCminc=0.2; % thermal State of Charge (minimum)
echc=1.00; edchc=1.00;
cchmax=1; % thermal max allowed charging power in kW
cdchmax=1; % thermal max allowed discharging power in kW
cselfdch=0.005; %self-discharging parameters
CScost=200; CSmc=0.3; CSyr=10; %TES unit cost, maintenance and life span
CSdgcost=0.1; %TES charging and discharging degradation cost
CSrep=(Ta/CSyr)-1;

%============== hygrogen storage variables==========
SOCminh=0.2; % Battery State of Charge (minimum)
echh2=0.98; edchh2=0.98;
h2chmax=1; % Battery max allowed charging power in kW
h2dchmax=1; % Battery max allowed discharging power in kW
selfdchh=0.001; %self-discharging parameters
HScost=200; HSmc=0.3; HSyr=10; %TES unit cost, maintenance and life span
HSdgcost=0.1; %TES charging and discharging degradation cost
HSrep=(Ta/HSyr)-1;

%=========environmental parameters======
CO2g=20; CO2grid=10; %carbon emission rate weight in kg/kw of NG and grid
cCO2=1.2; %unit penalty cost of carbon emission

%======initial time variables for 24hrs======
dt=1;       % Time step
nHours=numel(LD);
Time=(1:nHours)';
idxHr2Toend=2:nHours;

%initiate the problem
mesprob=optimproblem;

%initiate the variables

 % Wind Generation Constraint
 if WindVelocity(1:nHours)<Vci 
   Pwind=0;
  elseif WindVelocity(1:nHours)>Vco
      Pwind=0;
  elseif Vr<WindVelocity(1:nHours)<Vco
      Pwind=5;
  elseif Vci<WindVelocity(1:nHours)<Vr
      Pwind=5*((WindVelocity(1:nHours)-Vci)/(Vr-Vci));
 end

%============continuous variables=============
Pashp=optimvar('Pashp',nHours,1,'LowerBound',0);
Pecin=optimvar('Pecin',nHours,1,'LowerBound',0);
Pgshp=optimvar('Pgshp',nHours,1,'LowerBound',0);
gbGsh=optimvar('gbGsh',nHours,1,'LowerBound',0);
gbGas=optimvar('gbGas',nHours,1,'LowerBound',0);
Qcp=optimvar('Qcp',nHours,1,'LowerBound',0);
iceg=optimvar('iceg',nHours,1,'LowerBound',0);
grce=optimvar('grce',nHours,1,'LowerBound',0);
h1=optimvar('h1',nHours,1,'LowerBound',0);
h2=optimvar('h2',nHours,1,'LowerBound',0);
h3=optimvar('h3',nHours,1,'LowerBound',0);
Gpmofc=optimvar('Gpmofc',nHours,1,'LowerBound',0);
GgSOF=optimvar('GgSOF',nHours,1,'LowerBound',0);
Pelez=optimvar('Pelez',nHours,1,'LowerBound',0);
QAC=optimvar('QAC',nHours,1,'LowerBound',0);

%==========equipment capacity range===========
CPW=optimvar('CPW','LowerBound',5,'UpperBound',50);
CPV=optimvar('CPV','LowerBound',2.5,'UpperBound',50); % for PV
CPTS=optimvar('CPTS','LowerBound',5,'UpperBound',50);
CPRCE=optimvar('CPRCE','LowerBound',5,'UpperBound',30);
CPICE=optimvar('CPICE','LowerBound',5,'UpperBound',15);
CPEM=optimvar('CPEM','LowerBound',5,'UpperBound',15);
CPSOF=optimvar('CPSOF','LowerBound',5,'UpperBound',15);
CPMOF=optimvar('CPMOF','LowerBound',5,'UpperBound',15);
CQashp=optimvar('CQashp','LowerBound',5,'UpperBound',15);
CQwshp=optimvar('CQwshp','LowerBound',5,'UpperBound',15);
CPAC=optimvar('CPAC','LowerBound',5,'UpperBound',15);
CPEC=optimvar('CPEC','LowerBound',5,'UpperBound',15);
CPELEZ=optimvar('CPELEZ','LowerBound',5,'UpperBound',15);
CQgshp=optimvar('CQgshp','LowerBound',5,'UpperBound',15);
CQgb=optimvar('CQgb','LowerBound',5,'UpperBound',15);

Pgrid=optimvar('Pgrid',nHours,1,'LowerBound',0,'UpperBound',5); % for grid power
TSA=optimvar('TSA',nHours,1,'LowerBound',0,'UpperBound',300); % for grid power
Psell=optimvar('Psell',nHours,1,'LowerBound',0,'UpperBound',inf); % for selling power to grid
Et=optimvar('Et',nHours,1,'LowerBound',0,'UpperBound',10);
Eq=optimvar('Eq',nHours,1,'LowerBound',0,'UpperBound',10);
Ec=optimvar('Ec',nHours,1,'LowerBound',0,'UpperBound',10);
Eh2=optimvar('Eh2',nHours,1,'LowerBound',0,'UpperBound',10);
Pbch=optimvar('PbCh',nHours,1,'LowerBound',0,'UpperBound',1); % for battery charge
Pbdch=optimvar('Pbdch',nHours,1,'LowerBound',0,'UpperBound',1); % for battery discharge
qch=optimvar('qch',nHours,1,'LowerBound',0,'UpperBound',1); % for TES charge
qdch=optimvar('qdch',nHours,1,'LowerBound',0,'UpperBound',1); % for TES discharge
h2ch=optimvar('h2Ch',nHours,1,'LowerBound',0,'UpperBound',1); % for H2S charge
h2dch=optimvar('h2dch',nHours,1,'LowerBound',0,'UpperBound',1); % for H2S discharge
cch=optimvar('cch',nHours,1,'LowerBound',0,'UpperBound',1); % for cold  storage charge
cdch=optimvar('cdch',nHours,1,'LowerBound',0,'UpperBound',1); % for cold storage discharge
NOMb=optimvar('NOMb','LowerBound',5,'UpperBound',50); % nominal capacity of storage
NOMq=optimvar('NOMq','LowerBound',5,'UpperBound',50); % nominal capacity of storage
NOMc=optimvar('NOMc','LowerBound',5,'UpperBound',50); % nominal capacity of storage
NOMh=optimvar('NOMh','LowerBound',5,'UpperBound',50); % nominal capacity of storage


% ===================Initaiate the binary variables for selection==========
cp1=optimvar('cp1','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp2=optimvar('cp2','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp3=optimvar('cp3','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp4=optimvar('cp4','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp5=optimvar('cp5','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp6=optimvar('cp6','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp7=optimvar('cp7','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp8=optimvar('cp8','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cp9=optimvar('cp9','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq1=optimvar('cq1','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq2=optimvar('cq2','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq3=optimvar('cq3','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq4=optimvar('cq4','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
%========================operation binary variables=============
isONbch=optimvar('isONbch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge
isONbdch=optimvar('isONbdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge
isONqch=optimvar('isONqch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); 
isONqdch=optimvar('isONqdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
isONcch=optimvar('isONcch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); 
isONcdch=optimvar('isONcdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
isONh2ch=optimvar('isONh2Ch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); 
isONh2dch=optimvar('isONh2dch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
w1=optimvar('w1',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1); %binary to buy power
w2=optimvar('w2',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1); %binary to sell power
b1=optimvar('b1',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
b2=optimvar('b2',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
b3=optimvar('b3',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
b4=optimvar('b4',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
npv=optimvar('npv',nHours,1,'Type','integer','LowerBound',0,'UpperBound',5);  %number of pv to be installed
nwt=optimvar('nwt',nHours,1,'Type','integer','LowerBound',0,'UpperBound',5);    %number of wt to be installed

                %==============battery model============
mesprob.Constraints.Batstoredele=selfdchb*Et(idxHr2Toend-1,:)+(echb.*Pbch(idxHr2Toend,:)*dt)-...
    ((Pbdch(idxHr2Toend,:)*dt)./edchb)==Et(idxHr2Toend,:);
mesprob.Constraints.SOCb=Et(Time(1),:)<=Et(idxHr2Toend(end),:);
mesprob.Constraints.Etout=Et<= NOMb;
                %===============thermal storage model===========
mesprob.Constraints.thermalstored=(1-qselfdch)*Eq(idxHr2Toend-1,:)+(echq.*(qch(idxHr2Toend,:)*dt))-...
    ((qdch(idxHr2Toend,:)*dt)./edchq)==Eq(idxHr2Toend,:);
mesprob.Constraints.SOCq=Eq(Time(1),:)<=Eq(idxHr2Toend(end),:);
mesprob.Constraints.Eqout=Eq<= NOMq;
                %==========cold storage model==========
mesprob.Constraints.coldstored=(1-cselfdch)*Ec(idxHr2Toend-1,:)+(echc.*(cch(idxHr2Toend,:)*dt))...
    -((cdch(idxHr2Toend,:)*dt)./edchc)==Ec(idxHr2Toend,:);
mesprob.Constraints.SOCc=Ec(Time(1),:)<=Ec(idxHr2Toend(end),:);
mesprob.Constraints.Ecout=Ec<= NOMc;

%================hydrogen storage model====================================
mesprob.Constraints.hydeogenstored=selfdchh*Eh2(idxHr2Toend-1,:)+(echh2*(h2ch(idxHr2Toend,:)*dt))...
    -((h2dch(idxHr2Toend,:)*dt)./edchh2)==Eh2(idxHr2Toend,:);
mesprob.Constraints.SOCh=Eh2(Time(1),:)<=Eh2(idxHr2Toend(end),:);
mesprob.Constraints.Ehout=Eh2 <= NOMh;

%===================all equality constraints===============================
PVt=npv*PR*fR.*(SolarIrradiation'/GS);
PWt=nwt*Pwind;                              %wind power output
Pts=TSeff*TSA.*SolarIrradiation';           %thermal collector
Prce=gLHV*grce*RCEele;                      %electric output of RCE
Pice=(gLHV*iceelec*iceg);                   %electric output of ICE
Ppemfc=hLHV*h1*Pemeff;
Psofc=((hLHV*h2.*b1)+(gLHV*GgSOF.*b2))*sofeff;
Pmofc=((hLHV*h3.*b3)+(gLHV*Gpmofc.*b4))*Mofceff;
Helz=(elez/hLHV).*Pelez;                         %electrolyser output
AC=copAC*QAC;                                    %absorption chiller output
EC=copEC*Pecin;                                %electric chiller output
Qashp=Ashpeff*Pashp;                            %air source hp
Qgshp=gbheff*gbGsh*gLHV;                        %gas source hp
Qwshp=gshpef*Pgshp;                             %ground source hp
Qgb=gbeff*gbGas*gLHV;                           %gas boiler output

%=================capacity sizing==========================================
mesprob.Constraints.consh=Qgb<=CQgb; mesprob.Constraints.conso=Qwshp<=CQwshp; 
mesprob.Constraints.consa=PVt<=CPV; mesprob.Constraints.consi=PWt<=CPW; 
mesprob.Constraints.consb=Pts<=CPTS; mesprob.Constraints.consj=Prce<=CPRCE; 
mesprob.Constraints.consc=Pice<=CPICE; mesprob.Constraints.consk=Ppemfc<=CPEM; 
mesprob.Constraints.consd=Psofc<=CPSOF; mesprob.Constraints.consl=Pmofc<=CPMOF; 
mesprob.Constraints.conse=Helz<=CPELEZ; mesprob.Constraints.consm=AC<=CPAC; 
mesprob.Constraints.consf=EC<=CPEC; mesprob.Constraints.consn=Qashp<=CQashp; 
mesprob.Constraints.consg=Qgshp<=CQgshp;

Qpemfc=Ppemfc*((1-Pemeff-phloss)/Pemeff);
Qsofc=Psofc*((1-sofeff-sofhloss)/sofeff);
Qmofc=Pmofc*((1-Mofceff-mofhloss)/Mofceff);
Qrce=Prce*((1-RCEele-RCEloss)/RCEele);           %heat output of RCE
Qice=Pice*((1-iceelec-iceloss)/iceelec);         %heat output of ICE
Qfc=Qpemfc+Qsofc+Qmofc;                           %total heat from FC
Qchp=Qrce+Qice;                                    %total heat output from gas turbine
NG=Gpmofc+GgSOF+gbGas+gbGsh;                     %total natural gas demand

%====================all inequality constraints============================
mesprob.Constraints.econs1=Pbch<=PbChmax*isONbch; % Battery Limit of allowed Charging 
mesprob.Constraints.econs2=Pbdch<=PbDischmax*isONbdch; % Limit of allowed Discharging
mesprob.Constraints.econs3=isONbch+isONbdch<=1; % Forbidden the battery's charging and discharging simultaneously
mesprob.Constraints.econs4=qch<=qchmax*isONqch; % TES Limit of allowed Charging 
mesprob.Constraints.econs5=qdch<=qdchmax*isONqdch; % Limit of allowed Discharging
mesprob.Constraints.econs6=isONqch+isONqdch<=1;
mesprob.Constraints.ccons7=cch<=cchmax*isONcch; % CS Limit of allowed Charging 
mesprob.Constraints.econs8=cdch<=cdchmax*isONcdch; % CS Limit of allowed Discharging
mesprob.Constraints.econs9=isONcch+isONcdch<=1;
mesprob.Constraints.econs10=h2ch<=h2chmax*isONh2ch; % H2 storage Limit of allowed Charging 
mesprob.Constraints.econs11=h2dch<=h2dchmax*isONh2dch; % Limit of allowed Discharging
mesprob.Constraints.econs12=isONh2ch+isONh2dch<=1;
mesprob.Constraints.econs13=w1+w2<=1; %forbidden buy and selling same time.
mesprob.Constraints.Bsel=b1+b2<=1;
mesprob.Constraints.Bse2=b3+b4<=1;


%===============Carbon emission constraint==============
mesprob.Constraints.Cco2=sum((Gpmofc+GgSOF+gbGas+gbGsh)*dt*CO2g,1)+sum((Pgrid*dt)*CO2grid,1)<=...
    v.*((sum(LD'+HD'+CD')./Tloss).*1/grideff); %where v varies btw 0 and 1

%===============power balance=================

mesprob.Constraints.powerbalance=PVt(Time,:)+PWt(Time,:)+Prce(Time,:)+Pice(Time,:)...
    +Psofc(Time,:)+Pmofc(Time,:)+Ppemfc(Time,:)-Pbch(Time,:)+Pbdch(Time,:)+(w1.*Pgrid(Time,:))-(w2.*Psell(Time,:))==Pecin(Time,:)+Pashp(Time,:)+Pgshp(Time,:)+Pelez(Time,:)+LD';
%=============Heat balance==================
mesprob.Constraints.heatbalance= Qfc(Time,:)+Qchp(Time,:)+Qashp(Time,:)+Qwshp(Time,:)+Qgb(Time,:)+Qgshp(Time,:)+Pts(Time,:)-qch+qdch-HD'==0;
%==============hydrogen balance==============
mesprob.Constraints.hydrogenbalance=h1+h2+h3-Helz(Time,:)-h2ch+h2dch==0;
%==============cold balance==================
mesprob.Constraints.coldbalance=AC(Time,:)+EC(Time,:)-cch+cdch-CD'==0;

%=================objective function formulation================
%the total cost consist of annnualised investment cost, operation, maintenance and
%replacement cost% %annualised investment cost of each equipment
%==============investment cost of each equipment============
Cpvinv=pvcost*CPV*cp1; Cpwinv=pwcost*CPW*cp2; Cprceinv=rcecost*CPRCE*cp3; CpPTSinv=CPTS*TScost*cp9;
Cpiceinv=icecost*CPICE*cp4; Cpeminv=pemcost*CPEM*cp5; Cpsofinv=sofcost*CPSOF*cp6; Cpmofinv=mofcost*CPMOF*cp7;
Cpelzinv=elzcost*CPELEZ*cp8; CpBinv=BEcost*NOMb;
Cptesinv=TEScost*NOMq; Cpashinv=ashcost*CQashp*cq1; Cpwshpinv=wshpcost*CQwshp*cq2; Cpgbinv=gbcost*CQgb*cq3;
Cpgshpinv=gshpcost*CQgshp*cq4; CpACinv=CPAC*ACcost; CpECinv=CPEC*ECcost;
CpHSinv=HScost*NOMh; CpCSinv=NOMc*CScost;

Cinv=(Cpvinv+Cpwinv+Cprceinv+Cpiceinv+Cpeminv+Cpsofinv+Cpmofinv+Cpelzinv...
    +CpBinv+Cptesinv+Cpashinv+Cpwshpinv+Cpgbinv+Cpgshpinv+CpACinv+CpECinv+CpHSinv+CpCSinv+CpPTSinv)*AF;

%==============Maintenance and operational cost of equipment=============
Cm1=(Cprceinv*rcemc)+(Cpiceinv*icemc)+(Cpvinv*pvmc)+(Cpwinv*pwmc)+(Cpeminv*pemmc)+(Cpsofinv*sofmc)+(Cpmofinv*mofmc)...
    +(Cpelzinv+elzmc)+(CpPTSinv*TSmc)+(CpBinv*BEmc);
Cm2=(Cptesinv*TESmc)+(Cpashinv*ashmc)+(Cpwshpinv*wshpmc)+(Cpgbinv*gbmc)+(CpACinv*ACmc)+(CpECinv*ECmc)...
    +(CpHSinv*HSmc)+(CpCSinv*CSmc)+(Cpgshpinv*gshpmc);
   Cmain=Cm1+Cm2;

%===============Operation cost on fuel and gas==================
NGcost=sum(((NG*dt).*Cgas),1); gridcost=sum(((Pgrid*dt).*Cgrid'),1); sellcost=sum(((Psell*dt).*SellingCost),1);

%==============energy storages degradation cost=============
EESc1=sum(((Pbdch+Pbch)*dt)*BEdgcost,1)+sum(((qch+qdch)*dt)*TESdgcost,1);
EESc2=sum(((cch+cdch)*dt)*CSdgcost,1)+sum(((h2ch+h2dch)*dt)*HSdgcost,1);
Cop=(NGcost+gridcost+EESc1+EESc2-sellcost); %annual total operation cost

%===================replacement Cost========================
rep1=(Cprceinv*rcerep)+(Cpiceinv*icerep)+(Cpvinv*pvrep)+(Cpwinv*wtrep)+(Cpeminv*pemrep)+(Cpsofinv*sofrep)...
    +(Cpmofinv*mofrep)+(Cpelzinv*elzrep)+(CpBinv*Brep);
rep2=(Cptesinv*TESrep)+(Cpashinv*ashprep)+( Cpwshpinv*wshprep)+(Cpgbinv*gbrep)+(Cpgshpinv*gshphrep)...
    +(CpACinv*ACrep)+(CpECinv*ECrep)+(CpHSinv*HSrep)+(CpCSinv*CSrep)+(CpPTSinv*TSrep);
Creplacement=(rep1+rep2)*AF;

%===========set the objective function================
Totalcost=Cinv+Cmain+Creplacement+Cop;
mesprob.Objective=Totalcost;

 problem = prob2struct(mesprob);
 f = problem.objective;
Aeq = problem.Aeq;
beq =problem.beq;
Aineq = problem.Aineq;
bineq = problem.bineq;
lb = problem.lb;
ub = problem.ub;

% options for the optimization algorithm, here we set the max time it can run for
options = optimoptions('intlinprog','MaxTime',100000);
% call the optimization solver to find the best solution
[sol,TotalCost,exitflag,output] = solve(mesprob,'options',options);

[Matt J]

but the intlingprog identify the problem as non-linear, please how can i rectify it

Are you saying you don't believe the problem to be nonlinear?

[tobi alabi]

the problem is not non-linear, is just contain some binary variables

[tobi alabi]

@Matt J...the model is giving me this error:

Nonlinear problem with integer variables not supported.

please how i rectify it?

[Walter Roberson]

What would be some sample inputs?

[tobi alabi]

v=1; Ta=25.

Thanks Walter

[tobi alabi]

@walter @Matt J. I just modified the objective function to avoid the multiplication of two decision variables in the objective function, but i still have the same error. Please, how can correct this:

The modified code is;

%============= program for multi-energy system MES sizing=============

% this program consider optimal selection from the sub-categories of MES

% equipment, i.e CHP, PV, HP, GB. P2G, FUEL CELL

clc

% Load all the datas

LD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...

4 2.5 2.3 0.1]; % Home Load Data for 24 hrs in kW

HD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...

4 2.5 2.3 0.1]; % Heat Load Data for 24 hrs in kW

CD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...

4 2.5 2.3 0.1]; % cooling---------- Load Data for 24 hrs in kW

SolarIrradiation=[0 0 0 0 0.1 0.3 0.5 0.65 0.8 0.9 1 0.95 0.9 0.8 0.65...

0.4 0.3 0.08 0 0 0 0 0 0]; % for 24 hrs in kW/m^2

WindVelocity=repmat(5,1,24); % for 24 hrs in m/sec

Cgrid=[0.033 0.027 0.020 0.017 0.017 0.029 0.033 0.054 0.215...

0.572 0.572 0.572 0.215 0.572 0.286 0.279 0.086 0.059 0.050 0.061...

0.181 0.077 0.043 0.037]; % all prices and costs are in Euros

SellingCost=1.1; %feed in tarrif to the grid

Cgas=0.5;

% ==========PreDefined variables=================

%==========economy parameters

a=0.06; %interest rate, residual value and proposed lifecycle

v=input('please input the carbon emission target');

Ta=input('please input the projected service life of MES equipment');

Tloss=0.1; grideff=0.80;

AF1=a*((1+a)^Ta); af1=((1+a)^Ta)-1; %annuity factor

AF = AF1/af1;

%====== parameters variables for ICE and Reciprocating combustion engine

gLHV=25.3; %low heating value of natural gas in kw/m3

RCEele=0.75; iceelec=0.65; %electrical efficiency of the two CHP

RCEloss=0.03; iceloss=0.04; %heat loss factor of the two CHP

rceyr=15; iceyr=20; %life span of each CHP

rcecost=200; icecost=250; %unit cost of investment of each CHP

rcemc=0.3; icemc=0.4; % proportion of unit cost for maintenance

Prcemin=5; Prcemax=15; Picemin=5; Picemax=15;

rcerep=(Ta/rceyr)-1; icerep=(Ta/iceyr)-1;

%======= parameters variables for PEMFC, SOFC and MOFC===============

hLHV=33.3; %low heating value of hydrogen gas in kw/m3

Pemeff= 0.7; sofeff=0.8; Mofceff= 0.88; %electrical efficiency of FC

phloss=0.3; sofhloss=0.4; mofhloss=0.5; %heat loss of considered FC

Pemmin=5; Pemmax=15; Psomin=5; Psomax=15; Pmomin=5; Pmomax=15;

pemyr=15; sofyr=20; mofyr=25; %life span of each FC

pemcost=250; sofcost=300; mofcost=350; %unit cost of investment

pemmc=0.03; sofmc=0.04; mofmc=0.05; % proportion of unit cost for maintenance

pemrep=(Ta/pemyr)-1; sofrep=(Ta/sofyr)-1; mofrep=(Ta/mofyr)-1;

% ==============parameters variables for electrolyser============

elez=10; %specific electrical energy demand of electrolyser

eleyr=15; %life span of electrolyser

elzcost=300; %unit cost of elz investment

elzmc=0.03; %proportion of unit cost for maintenance

Hemax=15; Hemin=5;

elzrep=(Ta/eleyr)-1;

% =============photovoltaic parameters=============

PR=270; % rated power of selected pv module

fR=0.8; %derating factor of pv module

TSC=25; %standard temperature of pv cell

GS=1000; %rated solar radiation at kw/m2

pvyr=25; %pv life span of pv

pvcost=200; %pv unit invesmtent cost

pvmc=0.05; %proportion of investment cost for annual maintenance

PVmin=5; PVmax=10; pvrep=(Ta/pvyr)-1;

% ===========Solar Thermal collector parameters============

TSeff=0.65; %collector efficiency and installed area available

TSyr=25; TScost=200; TSmc=0.03; %collector cost, mcost and lifespan

TSAmin=0; TSAmax=300; TSrep=(Ta/TSyr)-1;

%============wind turbine parameters=============

Vci=4; % Cut-in Speed for wind turbine in m/sec

Vco=24; % Cut-out Speed for wind turbine in m/sec

Vr=14; % rated speed for wind turbine in m/sec

pw=5; %rated wind turbine power output

pwcost=250; %unit invetment cost of pv

pwmc=0.03; wtyr=20; %wind turbine mc cost and life span

wtrep=(Ta/wtyr)-1;

%===============parameters of heating equipment considered=================

Ashpeff=1.3; gshpef=1.4; %efficiency of air source & ground source hp

gbeff=0.92; gbheff=0.99; %efficiency of gas boiler and gas source hp

ashpyr=15; wshpyr=20; gbyr=15; gshpyr=15; %lifespan of each heat pump

ashcost=100; wshpcost=150; gbcost=200; gshpcost=250; %unit cost of hp

ashmc=0.03; wshpmc=0.04; gbmc=0.03; gshpmc=0.04; %maintenance cost

ashprep=(Ta/ashpyr)-1; gbrep=(Ta/gbyr)-1; wshprep=(Ta/wshpyr)-1;

gshphrep=(Ta/gshpyr)-1;

%===============parameters of cooling equipment considered=================

copAC=1.3; copEC=1.6; %coff. of performance of absorption and elec. chiller

ACyr=15; ECyr=20; %lifespan of chillers

ACcost=200; ECcost=250; %unit investment cost

ACmc=0.04; ECmc=0.03; %proportion of maintenance cost from investment cost

ACrep=(Ta/ACyr)-1; ECrep=(Ta/ECyr)-1;

% Energy storage variables

%======== electrical storage using electrochemical (battery)==============

echb=0.96; edchb=0.96; %charging abd discharching efficiency of Battery

PbChmax=1; % Battery max allowed charging power in kW

PbDischmax=1; % Battery max allowed discharging power in kW

selfdchb=0.001; %self-discharging parameters

BEcost=200; BEmc=0.03; BEyr=10; %BES unit cost, maintenance and life span

BEdgcost=0.1; %BES charging and discharging degradation cost

Brep=(Ta/BEyr)-1;

% ==========Thermal storage variables===========

qchmax=1; % thermal max allowed charging power in kW

qdchmax=1; % thermal max allowed discharging power in kW

qselfdch=0.005; %self-discharging parameters

echq=1.00; edchq=1.00; %charging abd discharching efficiency of Battery

tmin=10; tmax=35; %thermal storage temperature

TEScost=200; TESmc=0.3; TESyr=10; %TES unit cost, maintenance and life span

TESdgcost=0.1; %TES charging and discharging degradation cost

TESrep=(Ta/TESyr)-1;

% =============cold storage variables=============

SOCminc=0.2; % thermal State of Charge (minimum)

echc=1.00; edchc=1.00;

cchmax=1; % thermal max allowed charging power in kW

cdchmax=1; % thermal max allowed discharging power in kW

cselfdch=0.005; %self-discharging parameters

CScost=200; CSmc=0.3; CSyr=10; %TES unit cost, maintenance and life span

CSdgcost=0.1; %TES charging and discharging degradation cost

CSrep=(Ta/CSyr)-1;

%============== hygrogen storage variables==========

SOCminh=0.2; % Battery State of Charge (minimum)

echh2=0.98; edchh2=0.98;

h2chmax=1; % Battery max allowed charging power in kW

h2dchmax=1; % Battery max allowed discharging power in kW

selfdchh=0.001; %self-discharging parameters

HScost=200; HSmc=0.3; HSyr=10; %TES unit cost, maintenance and life span

HSdgcost=0.1; %TES charging and discharging degradation cost

HSrep=(Ta/HSyr)-1;

%=========environmental parameters======

CO2g=20; CO2grid=10; %carbon emission rate weight in kg/kw of NG and grid

cCO2=1.2; %unit penalty cost of carbon emission

%======initial time variables for 24hrs======

dt=1; % Time step

nHours=numel(LD);

Time=(1:nHours)';

idxHr2Toend=2:nHours;

%initiate the problem

mesprob=optimproblem;

%initiate the variables

% Wind Generation Constraint

if WindVelocity(1:nHours)<Vci

Pwind=0;

elseif WindVelocity(1:nHours)>Vco

Pwind=0;

elseif Vr<WindVelocity(1:nHours)<Vco

Pwind=5;

elseif Vci<WindVelocity(1:nHours)<Vr

Pwind=5*((WindVelocity(1:nHours)-Vci)/(Vr-Vci));

end

%========equipment min-max limits==========

CPWmin=5; CPWmax=50; CPVmin=5; CPVmax=50;

CPTSmin=5; CPTSmax=50; CPRCEmin=5; CPRCEmax=50;

CPICEmin=5; CPICEmax=15; CPEMmin=5; CPEMmax=15;

CPSOFmin=5; CPSOFmax=15; CPMOFmin=5; CPMOFmax=15;

CQashpmin=5; CQashpmax=15; CQwshpmin=5; CQwshpmax=15;

CPACmin=5; CPACmax=15; CPECmin=5; CPECmax=15;

CQgshpmin=5; CQgshpmax=15; CQgbmin=5; CQgbmax=15;

CPELEZmin=5; CPELEZmax=15;

%============continuous variables=============

Pashp=optimvar('Pashp',nHours,1,'LowerBound',0);

Pecin=optimvar('Pecin',nHours,1,'LowerBound',0);

Pgshp=optimvar('Pgshp',nHours,1,'LowerBound',0);

gbGsh=optimvar('gbGsh',nHours,1,'LowerBound',0);

gbGas=optimvar('gbGas',nHours,1,'LowerBound',0);

Qcp=optimvar('Qcp',nHours,1,'LowerBound',0);

iceg=optimvar('iceg',nHours,1,'LowerBound',0);

grce=optimvar('grce',nHours,1,'LowerBound',0);

h1=optimvar('h1',nHours,1,'LowerBound',0);

h2=optimvar('h2',nHours,1,'LowerBound',0);

h3=optimvar('h3',nHours,1,'LowerBound',0);

Gpmofc=optimvar('Gpmofc',nHours,1,'LowerBound',0);

GgSOF=optimvar('GgSOF',nHours,1,'LowerBound',0);

Pelez=optimvar('Pelez',nHours,1,'LowerBound',0);

QAC=optimvar('QAC',nHours,1,'LowerBound',0);

%==========equipment capacity range===========

CPW=optimvar('CPW','LowerBound',5,'UpperBound',50);

CPV=optimvar('CPV','LowerBound',2.5,'UpperBound',50); % for PV

CPTS=optimvar('CPTS','LowerBound',5,'UpperBound',50);

CPRCE=optimvar('CPRCE','LowerBound',5,'UpperBound',30);

CPICE=optimvar('CPICE','LowerBound',5,'UpperBound',15);

CPEM=optimvar('CPEM','LowerBound',5,'UpperBound',15);

CPSOF=optimvar('CPSOF','LowerBound',5,'UpperBound',15);

CPMOF=optimvar('CPMOF','LowerBound',5,'UpperBound',15);

CQashp=optimvar('CQashp','LowerBound',5,'UpperBound',15);

CQwshp=optimvar('CQwshp','LowerBound',5,'UpperBound',15);

CPAC=optimvar('CPAC','LowerBound',5,'UpperBound',15);

CPEC=optimvar('CPEC','LowerBound',5,'UpperBound',15);

CPELEZ=optimvar('CPELEZ','LowerBound',5,'UpperBound',15);

CQgshp=optimvar('CQgshp','LowerBound',5,'UpperBound',15);

CQgb=optimvar('CQgb','LowerBound',5,'UpperBound',15);

Pgrid=optimvar('Pgrid',nHours,1,'LowerBound',0,'UpperBound',5); % for grid power

TSA=optimvar('TSA',nHours,1,'LowerBound',0,'UpperBound',300); % for grid power

Psell=optimvar('Psell',nHours,1,'LowerBound',0,'UpperBound',inf); % for selling power to grid

Et=optimvar('Et',nHours,1,'LowerBound',0,'UpperBound',10);

Eq=optimvar('Eq',nHours,1,'LowerBound',0,'UpperBound',10);

Ec=optimvar('Ec',nHours,1,'LowerBound',0,'UpperBound',10);

Eh2=optimvar('Eh2',nHours,1,'LowerBound',0,'UpperBound',10);

Pbch=optimvar('PbCh',nHours,1,'LowerBound',0,'UpperBound',1); % for battery charge

Pbdch=optimvar('Pbdch',nHours,1,'LowerBound',0,'UpperBound',1); % for battery discharge

qch=optimvar('qch',nHours,1,'LowerBound',0,'UpperBound',1); % for TES charge

qdch=optimvar('qdch',nHours,1,'LowerBound',0,'UpperBound',1); % for TES discharge

h2ch=optimvar('h2Ch',nHours,1,'LowerBound',0,'UpperBound',1); % for H2S charge

h2dch=optimvar('h2dch',nHours,1,'LowerBound',0,'UpperBound',1); % for H2S discharge

cch=optimvar('cch',nHours,1,'LowerBound',0,'UpperBound',1); % for cold storage charge

cdch=optimvar('cdch',nHours,1,'LowerBound',0,'UpperBound',1); % for cold storage discharge

NOMb=optimvar('NOMb','LowerBound',5,'UpperBound',50); % nominal capacity of storage

NOMq=optimvar('NOMq','LowerBound',5,'UpperBound',50); % nominal capacity of storage

NOMc=optimvar('NOMc','LowerBound',5,'UpperBound',50); % nominal capacity of storage

NOMh=optimvar('NOMh','LowerBound',5,'UpperBound',50); % nominal capacity of storage

% ===================Initaiate the binary variables for selection==========

cp1=optimvar('cp1','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq

cp2=optimvar('cp2','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq

cp3=optimvar('cp3','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq

cp4=optimvar('cp4','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq

cp5=optimvar('cp5','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq

cp6=optimvar('cp6','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq

cp7=optimvar('cp7','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq

cp8=optimvar('cp8','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq

cp9=optimvar('cp9','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq

cq1=optimvar('cq1','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq

cq2=optimvar('cq2','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq

cq3=optimvar('cq3','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq

cq4=optimvar('cq4','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq

%========================operation binary variables=============

isONbch=optimvar('isONbch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge

isONbdch=optimvar('isONbdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge

isONqch=optimvar('isONqch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);

isONqdch=optimvar('isONqdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);

isONcch=optimvar('isONcch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);

isONcdch=optimvar('isONcdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);

isONh2ch=optimvar('isONh2Ch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);

isONh2dch=optimvar('isONh2dch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);

w1=optimvar('w1',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1); %binary to buy power

w2=optimvar('w2',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1); %binary to sell power

b1=optimvar('b1',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);

b2=optimvar('b2',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);

b3=optimvar('b3',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);

b4=optimvar('b4',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);

npv=optimvar('npv',nHours,1,'Type','integer','LowerBound',0,'UpperBound',5); %number of pv to be installed

nwt=optimvar('nwt',nHours,1,'Type','integer','LowerBound',0,'UpperBound',5); %number of wt to be installed

%==============battery model============

mesprob.Constraints.Batstoredele=selfdchb*Et(idxHr2Toend-1,:)+(echb.*Pbch(idxHr2Toend,:)*dt)-...

((Pbdch(idxHr2Toend,:)*dt)./edchb)==Et(idxHr2Toend,:);

mesprob.Constraints.SOCb=Et(Time(1),:)<=Et(idxHr2Toend(end),:);

mesprob.Constraints.Etout=Et<= NOMb;

%===============thermal storage model===========

mesprob.Constraints.thermalstored=(1-qselfdch)*Eq(idxHr2Toend-1,:)+(echq.*(qch(idxHr2Toend,:)*dt))-...

((qdch(idxHr2Toend,:)*dt)./edchq)==Eq(idxHr2Toend,:);

mesprob.Constraints.SOCq=Eq(Time(1),:)<=Eq(idxHr2Toend(end),:);

mesprob.Constraints.Eqout=Eq<= NOMq;

%==========cold storage model==========

mesprob.Constraints.coldstored=(1-cselfdch)*Ec(idxHr2Toend-1,:)+(echc.*(cch(idxHr2Toend,:)*dt))...

-((cdch(idxHr2Toend,:)*dt)./edchc)==Ec(idxHr2Toend,:);

mesprob.Constraints.SOCc=Ec(Time(1),:)<=Ec(idxHr2Toend(end),:);

mesprob.Constraints.Ecout=Ec<= NOMc;

%================hydrogen storage model====================================

mesprob.Constraints.hydeogenstored=selfdchh*Eh2(idxHr2Toend-1,:)+(echh2*(h2ch(idxHr2Toend,:)*dt))...

-((h2dch(idxHr2Toend,:)*dt)./edchh2)==Eh2(idxHr2Toend,:);

mesprob.Constraints.SOCh=Eh2(Time(1),:)<=Eh2(idxHr2Toend(end),:);

mesprob.Constraints.Ehout=Eh2 <= NOMh;

%===================all equality constraints===============================

PVt=npv*PR*fR.*(SolarIrradiation'/GS);

PWt=nwt*Pwind; %wind power output

Pts=TSeff*TSA.*SolarIrradiation'; %thermal collector

Prce=gLHV*grce*RCEele; %electric output of RCE

Pice=(gLHV*iceelec*iceg); %electric output of ICE

Ppemfc=hLHV*h1*Pemeff;

Psofc=((hLHV*h2.*b1)+(gLHV*GgSOF.*b2))*sofeff;

Pmofc=((hLHV*h3.*b3)+(gLHV*Gpmofc.*b4))*Mofceff;

Helz=(elez/hLHV).*Pelez; %electrolyser output

AC=copAC*QAC; %absorption chiller output

EC=copEC*Pecin; %electric chiller output

Qashp=Ashpeff*Pashp; %air source hp

Qgshp=gbheff*gbGsh*gLHV; %gas source hp

Qwshp=gshpef*Pgshp; %ground source hp

Qgb=gbeff*gbGas*gLHV; %gas boiler output

%=================capacity sizing==========================================

mesprob.Constraints.consa2=PVt<=CPV;

mesprob.Constraints.consa1=CPVmin*cp1<=CPV; mesprob.Constraints.consa=CPV<=CPVmax*cp1;

mesprob.Constraints.consi2=PWt<=CPW;

mesprob.Constraints.consi1=CPWmin*cp2<=CPW; mesprob.Constraints.consi=CPW<=CPWmax*cp2;

mesprob.Constraints.consb2=Pts<=CPTS;

mesprob.Constraints.consb1=CPTSmin*cp3<=CPTS; mesprob.Constraints.consb=CPTS<=CPTSmax*cp3;

mesprob.Constraints.consj2=Prce<=CPRCE;

mesprob.Constraints.consj1=CPRCEmin*cp4<=CPRCE; mesprob.Constraints.consj=CPRCE<=CPRCEmax*cp4;

mesprob.Constraints.consc2=Pice<=CPICE;

mesprob.Constraints.consc1=CPICEmin*cp5<=CPICE; mesprob.Constraints.consc=CPICE<=CPICEmax*cp5;

mesprob.Constraints.consk2=Ppemfc<=CPEM;

mesprob.Constraints.consk1=CPEMmin*cp6<=CPEM; mesprob.Constraints.consk=CPEM<=CPEMmax*cp6;

mesprob.Constraints.consd2=Psofc<=CPSOF;

mesprob.Constraints.consd1=CPSOFmin*cp7<=CPSOF; mesprob.Constraints.consd=CPSOF<=CPSOFmax*cp7;

mesprob.Constraints.consl=Pmofc<=CPMOF;

mesprob.Constraints.consl1=CPMOFmin*cp8<=CPMOF; mesprob.Constraints.consl=CPMOF<=CPMOFmax*cp8;

mesprob.Constraints.conse2=Helz<=CPELEZ;

mesprob.Constraints.conse1=CPELEZmin*cp9<=CPELEZ; mesprob.Constraints.conse=CPELEZ<=CPELEZmax*cp9;

mesprob.Constraints.consn2=Qashp<=CQashp;

mesprob.Constraints.consn1=CQashp<=CQashpmin*cq1; mesprob.Constraints.consn=CQashp<=CQashpmax*cq1;

mesprob.Constraints.consg2=Qgshp<=CQgshp;

mesprob.Constraints.consg1=CQgshpmin*cq2<=CQgshp; mesprob.Constraints.consg=CQgshp<=CQgshpmax*cq2;

mesprob.Constraints.consh2=Qgb<=CQgb;

mesprob.Constraints.consh1=CQgbmin*cq3<=CQgb; mesprob.Constraints.consh=CQgb<=CQgbmax*cq3;

mesprob.Constraints.conso2=Qwshp<=CQwshp;

mesprob.Constraints.conso1=CQwshpmin*cq4<=CQwshp; mesprob.Constraints.conso=CQwshp<=CQwshpmax*cq4;

mesprob.Constraints.consm=AC<=CPAC;

mesprob.Constraints.consf=EC<=CPEC;

Qpemfc=Ppemfc*((1-Pemeff-phloss)/Pemeff);

Qsofc=Psofc*((1-sofeff-sofhloss)/sofeff);

Qmofc=Pmofc*((1-Mofceff-mofhloss)/Mofceff);

Qrce=Prce*((1-RCEele-RCEloss)/RCEele); %heat output of RCE

Qice=Pice*((1-iceelec-iceloss)/iceelec); %heat output of ICE

Qfc=Qpemfc+Qsofc+Qmofc; %total heat from FC

Qchp=Qrce+Qice; %total heat output from gas turbine

NG=Gpmofc+GgSOF+gbGas+gbGsh; %total natural gas demand

%====================all inequality constraints============================

mesprob.Constraints.econs1=Pbch<=PbChmax*isONbch; % Battery Limit of allowed Charging

mesprob.Constraints.econs2=Pbdch<=PbDischmax*isONbdch; % Limit of allowed Discharging

mesprob.Constraints.econs3=isONbch+isONbdch<=1; % Forbidden the battery's charging and discharging simultaneously

mesprob.Constraints.econs4=qch<=qchmax*isONqch; % TES Limit of allowed Charging

mesprob.Constraints.econs5=qdch<=qdchmax*isONqdch; % Limit of allowed Discharging

mesprob.Constraints.econs6=isONqch+isONqdch<=1;

mesprob.Constraints.ccons7=cch<=cchmax*isONcch; % CS Limit of allowed Charging

mesprob.Constraints.econs8=cdch<=cdchmax*isONcdch; % CS Limit of allowed Discharging

mesprob.Constraints.econs9=isONcch+isONcdch<=1;

mesprob.Constraints.econs10=h2ch<=h2chmax*isONh2ch; % H2 storage Limit of allowed Charging

mesprob.Constraints.econs11=h2dch<=h2dchmax*isONh2dch; % Limit of allowed Discharging

mesprob.Constraints.econs12=isONh2ch+isONh2dch<=1;

mesprob.Constraints.econs13=w1+w2<=1; %forbidden buy and selling same time.

mesprob.Constraints.Bsel=b1+b2<=1;

mesprob.Constraints.Bse2=b3+b4<=1;

%===============Carbon emission constraint==============

mesprob.Constraints.Cco2=sum((Gpmofc+GgSOF+gbGas+gbGsh)*dt*CO2g,1)+sum((Pgrid*dt)*CO2grid,1)<=...

v.*((sum(LD'+HD'+CD')./Tloss).*1/grideff); %where v varies btw 0 and 1

%===============power balance=================

mesprob.Constraints.powerbalance=PVt(Time,:)+PWt(Time,:)+Prce(Time,:)+Pice(Time,:)...

+Psofc(Time,:)+Pmofc(Time,:)+Ppemfc(Time,:)-Pbch(Time,:)+Pbdch(Time,:)+(w1.*Pgrid(Time,:))-(w2.*Psell(Time,:))==Pecin(Time,:)+Pashp(Time,:)+Pgshp(Time,:)+Pelez(Time,:)+LD';

%=============Heat balance==================

mesprob.Constraints.heatbalance= Qfc(Time,:)+Qchp(Time,:)+Qashp(Time,:)+Qwshp(Time,:)+Qgb(Time,:)+Qgshp(Time,:)+Pts(Time,:)-qch+qdch-HD'==0;

%==============hydrogen balance==============

mesprob.Constraints.hydrogenbalance=h1+h2+h3-Helz(Time,:)-h2ch+h2dch==0;

%==============cold balance==================

mesprob.Constraints.coldbalance=AC(Time,:)+EC(Time,:)-cch+cdch-CD'==0;

%=================objective function formulation================

%the total cost consist of annnualised investment cost, operation, maintenance and

%replacement cost% %annualised investment cost of each equipment

%==============investment cost of each equipment============

Cpvinv=pvcost*CPV; Cpwinv=pwcost*CPW; Cprceinv=rcecost*CPRCE; CpPTSinv=CPTS*TScost;

Cpiceinv=icecost*CPICE; Cpeminv=pemcost*CPEM; Cpsofinv=sofcost*CPSOF; Cpmofinv=mofcost*CPMOF;

Cpelzinv=elzcost*CPELEZ; CpBinv=BEcost*NOMb;

Cptesinv=TEScost*NOMq; Cpashinv=ashcost*CQashp; Cpwshpinv=wshpcost*CQwshp; Cpgbinv=gbcost*CQgb;

Cpgshpinv=gshpcost*CQgshp; CpACinv=CPAC*ACcost; CpECinv=CPEC*ECcost;

CpHSinv=HScost*NOMh; CpCSinv=NOMc*CScost;

Cinv=(Cpvinv+Cpwinv+Cprceinv+Cpiceinv+Cpeminv+Cpsofinv+Cpmofinv+Cpelzinv...

+CpBinv+Cptesinv+Cpashinv+Cpwshpinv+Cpgbinv+Cpgshpinv+CpACinv+CpECinv+CpHSinv+CpCSinv+CpPTSinv)*AF;

%==============Maintenance and operational cost of equipment=============

Cm1=(Cprceinv*rcemc)+(Cpiceinv*icemc)+(Cpvinv*pvmc)+(Cpwinv*pwmc)+(Cpeminv*pemmc)+(Cpsofinv*sofmc)+(Cpmofinv*mofmc)...

+(Cpelzinv+elzmc)+(CpPTSinv*TSmc)+(CpBinv*BEmc);

Cm2=(Cptesinv*TESmc)+(Cpashinv*ashmc)+(Cpwshpinv*wshpmc)+(Cpgbinv*gbmc)+(CpACinv*ACmc)+(CpECinv*ECmc)...

+(CpHSinv*HSmc)+(CpCSinv*CSmc)+(Cpgshpinv*gshpmc);

Cmain=Cm1+Cm2;

%===============Operation cost on fuel and gas==================

NGcost=sum(((NG*dt).*Cgas),1); gridcost=sum(((Pgrid*dt).*Cgrid'),1); sellcost=sum(((Psell*dt).*SellingCost),1);

%==============energy storages degradation cost=============

EESc1=sum(((Pbdch+Pbch)*dt)*BEdgcost,1)+sum(((qch+qdch)*dt)*TESdgcost,1);

EESc2=sum(((cch+cdch)*dt)*CSdgcost,1)+sum(((h2ch+h2dch)*dt)*HSdgcost,1);

Cop=(NGcost+gridcost+EESc1+EESc2-sellcost); %annual total operation cost

%===================replacement Cost========================

rep1=(Cprceinv*rcerep)+(Cpiceinv*icerep)+(Cpvinv*pvrep)+(Cpwinv*wtrep)+(Cpeminv*pemrep)+(Cpsofinv*sofrep)...

+(Cpmofinv*mofrep)+(Cpelzinv*elzrep)+(CpBinv*Brep);

rep2=(Cptesinv*TESrep)+(Cpashinv*ashprep)+( Cpwshpinv*wshprep)+(Cpgbinv*gbrep)+(Cpgshpinv*gshphrep)...

+(CpACinv*ACrep)+(CpECinv*ECrep)+(CpHSinv*HSrep)+(CpCSinv*CSrep)+(CpPTSinv*TSrep);

Creplacement=(rep1+rep2)*AF;

%===========set the objective function================

Totalcost=Cinv+Cmain+Creplacement+Cop;

mesprob.Objective=Totalcost;

%===========options for the optimization algorithm====================

%options = cplexoptimset('cplex');

% options.display = 'on';

% ctype = 'CCCI';

%==========call the optimization solver to find the best solution=======

% options for the optimization algorithm, here we set the max time it can run for

options = optimoptions('intlinprog','MaxTime',100000);

% call the optimization solver to find the best solution

[sol,TotalCost,exitflag,output] = solve(mesprob,'options',options);

%showproblem(mesprob);

[Matt J]

Please put the code in a code formating block so that it is easier for us to copy/paste and test.

[tobi alabi]

%============= program for multi-energy system MES sizing=============
% this program consider optimal selection from the sub-categories of MES
% equipment, i.e CHP, PV, HP, GB. P2G, FUEL CELL
clc
% Load all the datas
LD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...
    4 2.5 2.3 0.1]; % Home Load Data for 24 hrs in kW
HD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...
    4 2.5 2.3 0.1]; % Heat Load Data for 24 hrs in kW
CD=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...
    4 2.5 2.3 0.1]; % cooling---------- Load Data for 24 hrs in kW
SolarIrradiation=[0 0 0 0 0.1 0.3 0.5 0.65 0.8 0.9 1 0.95 0.9 0.8 0.65...
    0.4 0.3 0.08 0 0 0 0 0 0]; % for 24 hrs in kW/m^2
WindVelocity=repmat(5,1,24); % for 24 hrs in m/sec
Cgrid=[0.033 0.027 0.020 0.017 0.017 0.029 0.033 0.054 0.215...
    0.572 0.572 0.572 0.215 0.572 0.286 0.279 0.086 0.059 0.050 0.061...
    0.181 0.077 0.043 0.037]; % all prices and costs are in Euros
SellingCost=1.1; %feed in tarrif to the grid
Cgas=0.5;
% ==========PreDefined variables=================
%==========economy parameters
a=0.06;  %interest rate, residual value and proposed lifecycle
v=input('please input the carbon emission target');
Ta=input('please input the projected service life of MES equipment');
Tloss=0.1; grideff=0.80;
AF1=a*((1+a)^Ta); af1=((1+a)^Ta)-1; %annuity factor
AF = AF1/af1;
%====== parameters variables for ICE and Reciprocating combustion engine
gLHV=25.3; %low heating value of natural gas in kw/m3
RCEele=0.75; iceelec=0.65; %electrical efficiency of the two CHP
RCEloss=0.03; iceloss=0.04; %heat loss factor of the two CHP
rceyr=15; iceyr=20; %life span of each CHP
rcecost=200; icecost=250; %unit cost of investment of each CHP
rcemc=0.3; icemc=0.4; % proportion of unit cost for maintenance
Prcemin=5; Prcemax=15; Picemin=5; Picemax=15;
rcerep=(Ta/rceyr)-1; icerep=(Ta/iceyr)-1;

%======= parameters variables for PEMFC, SOFC and MOFC===============
hLHV=33.3; %low heating value of hydrogen gas in kw/m3
Pemeff= 0.7; sofeff=0.8; Mofceff= 0.88; %electrical efficiency of FC
phloss=0.3; sofhloss=0.4; mofhloss=0.5; %heat loss of considered FC
Pemmin=5; Pemmax=15; Psomin=5; Psomax=15; Pmomin=5; Pmomax=15;
pemyr=15; sofyr=20; mofyr=25; %life span of each FC
pemcost=250; sofcost=300; mofcost=350; %unit cost of investment
pemmc=0.03; sofmc=0.04; mofmc=0.05; % proportion of unit cost for maintenance
pemrep=(Ta/pemyr)-1; sofrep=(Ta/sofyr)-1; mofrep=(Ta/mofyr)-1;

% ==============parameters variables for electrolyser============
elez=10; %specific electrical energy demand of electrolyser
eleyr=15; %life span of electrolyser
elzcost=300; %unit cost of elz investment
elzmc=0.03; %proportion of unit cost for maintenance
Hemax=15; Hemin=5;
elzrep=(Ta/eleyr)-1;

% =============photovoltaic parameters=============
PR=270; % rated power of selected pv module
fR=0.8; %derating factor of pv module
TSC=25; %standard temperature of pv cell
GS=1000; %rated solar radiation at kw/m2
pvyr=25; %pv life span of pv
pvcost=200; %pv unit invesmtent cost
pvmc=0.05; %proportion of investment cost for annual maintenance
PVmin=5; PVmax=10; pvrep=(Ta/pvyr)-1;

% ===========Solar Thermal collector parameters============
TSeff=0.65;  %collector efficiency and installed area available
TSyr=25; TScost=200; TSmc=0.03; %collector cost, mcost and lifespan
TSAmin=0; TSAmax=300; TSrep=(Ta/TSyr)-1;

%============wind turbine parameters=============
Vci=4; % Cut-in Speed for wind turbine in m/sec
Vco=24; % Cut-out Speed for wind turbine in m/sec
Vr=14; % rated speed for wind turbine in m/sec
pw=5; %rated wind turbine power output
pwcost=250; %unit invetment cost of pv
pwmc=0.03; wtyr=20; %wind turbine mc cost and life span
wtrep=(Ta/wtyr)-1;

%===============parameters of heating equipment considered=================
Ashpeff=1.3; gshpef=1.4; %efficiency of air source & ground source hp
gbeff=0.92; gbheff=0.99; %efficiency of gas boiler and gas source hp
ashpyr=15; wshpyr=20; gbyr=15; gshpyr=15; %lifespan of each heat pump
ashcost=100; wshpcost=150; gbcost=200; gshpcost=250; %unit cost of hp
ashmc=0.03; wshpmc=0.04; gbmc=0.03; gshpmc=0.04; %maintenance cost
ashprep=(Ta/ashpyr)-1; gbrep=(Ta/gbyr)-1; wshprep=(Ta/wshpyr)-1;
gshphrep=(Ta/gshpyr)-1;

%===============parameters of cooling equipment considered=================
copAC=1.3; copEC=1.6; %coff. of performance of absorption and elec. chiller
ACyr=15; ECyr=20; %lifespan of chillers
ACcost=200; ECcost=250; %unit investment cost
ACmc=0.04; ECmc=0.03; %proportion of maintenance cost from investment cost 
ACrep=(Ta/ACyr)-1; ECrep=(Ta/ECyr)-1;

            % Energy storage variables
%======== electrical storage using electrochemical (battery)==============
echb=0.96; edchb=0.96;       %charging abd discharching efficiency of Battery
PbChmax=1;                  % Battery max allowed charging power in kW
PbDischmax=1;                % Battery max allowed discharging power in kW
selfdchb=0.001;                 %self-discharging parameters
BEcost=200; BEmc=0.03; BEyr=10; %BES unit cost, maintenance and life span
BEdgcost=0.1;                   %BES charging and discharging degradation cost
Brep=(Ta/BEyr)-1;

% ==========Thermal storage variables===========
qchmax=1;                   % thermal max allowed charging power in kW
qdchmax=1;                   % thermal max allowed discharging power in kW
qselfdch=0.005;             %self-discharging parameters
echq=1.00; edchq=1.00;      %charging abd discharching efficiency of Battery
tmin=10; tmax=35;            %thermal storage temperature
TEScost=200; TESmc=0.3; TESyr=10; %TES unit cost, maintenance and life span
TESdgcost=0.1;               %TES charging and discharging degradation cost
TESrep=(Ta/TESyr)-1;
  
% =============cold storage variables=============
SOCminc=0.2; % thermal State of Charge (minimum)
echc=1.00; edchc=1.00;
cchmax=1; % thermal max allowed charging power in kW
cdchmax=1; % thermal max allowed discharging power in kW
cselfdch=0.005; %self-discharging parameters
CScost=200; CSmc=0.3; CSyr=10; %TES unit cost, maintenance and life span
CSdgcost=0.1; %TES charging and discharging degradation cost
CSrep=(Ta/CSyr)-1;

%============== hygrogen storage variables==========
SOCminh=0.2; % Battery State of Charge (minimum)
echh2=0.98; edchh2=0.98;
h2chmax=1; % Battery max allowed charging power in kW
h2dchmax=1; % Battery max allowed discharging power in kW
selfdchh=0.001; %self-discharging parameters
HScost=200; HSmc=0.3; HSyr=10; %TES unit cost, maintenance and life span
HSdgcost=0.1; %TES charging and discharging degradation cost
HSrep=(Ta/HSyr)-1;

%=========environmental parameters======
CO2g=20; CO2grid=10; %carbon emission rate weight in kg/kw of NG and grid
cCO2=1.2; %unit penalty cost of carbon emission

%======initial time variables for 24hrs======
dt=1;       % Time step
nHours=numel(LD);
Time=(1:nHours)';
idxHr2Toend=2:nHours;

%initiate the problem
mesprob=optimproblem;

%initiate the variables

 % Wind Generation Constraint
 if WindVelocity(1:nHours)<Vci 
   Pwind=0;
  elseif WindVelocity(1:nHours)>Vco
      Pwind=0;
  elseif Vr<WindVelocity(1:nHours)<Vco
      Pwind=5;
  elseif Vci<WindVelocity(1:nHours)<Vr
      Pwind=5*((WindVelocity(1:nHours)-Vci)/(Vr-Vci));
 end
 %========equipment min-max limits==========
 CPWmin=5; CPWmax=50; CPVmin=5; CPVmax=50;
CPTSmin=5; CPTSmax=50; CPRCEmin=5; CPRCEmax=50;
CPICEmin=5; CPICEmax=15; CPEMmin=5; CPEMmax=15;
CPSOFmin=5; CPSOFmax=15; CPMOFmin=5; CPMOFmax=15;
CQashpmin=5; CQashpmax=15; CQwshpmin=5; CQwshpmax=15;
CPACmin=5; CPACmax=15; CPECmin=5; CPECmax=15;
CQgshpmin=5; CQgshpmax=15; CQgbmin=5; CQgbmax=15;
CPELEZmin=5; CPELEZmax=15;
%============continuous variables=============
Pashp=optimvar('Pashp',nHours,1,'LowerBound',0);
Pecin=optimvar('Pecin',nHours,1,'LowerBound',0);
Pgshp=optimvar('Pgshp',nHours,1,'LowerBound',0);
gbGsh=optimvar('gbGsh',nHours,1,'LowerBound',0);
gbGas=optimvar('gbGas',nHours,1,'LowerBound',0);
Qcp=optimvar('Qcp',nHours,1,'LowerBound',0);
iceg=optimvar('iceg',nHours,1,'LowerBound',0);
grce=optimvar('grce',nHours,1,'LowerBound',0);
h1=optimvar('h1',nHours,1,'LowerBound',0);
h2=optimvar('h2',nHours,1,'LowerBound',0);
h3=optimvar('h3',nHours,1,'LowerBound',0);
Gpmofc=optimvar('Gpmofc',nHours,1,'LowerBound',0);
GgSOF=optimvar('GgSOF',nHours,1,'LowerBound',0);
Pelez=optimvar('Pelez',nHours,1,'LowerBound',0);
QAC=optimvar('QAC',nHours,1,'LowerBound',0);

%==========equipment capacity range===========
CPW=optimvar('CPW','LowerBound',5,'UpperBound',50);
CPV=optimvar('CPV','LowerBound',2.5,'UpperBound',50); % for PV
CPTS=optimvar('CPTS','LowerBound',5,'UpperBound',50);
CPRCE=optimvar('CPRCE','LowerBound',5,'UpperBound',30);
CPICE=optimvar('CPICE','LowerBound',5,'UpperBound',15);
CPEM=optimvar('CPEM','LowerBound',5,'UpperBound',15);
CPSOF=optimvar('CPSOF','LowerBound',5,'UpperBound',15);
CPMOF=optimvar('CPMOF','LowerBound',5,'UpperBound',15);
CQashp=optimvar('CQashp','LowerBound',5,'UpperBound',15);
CQwshp=optimvar('CQwshp','LowerBound',5,'UpperBound',15);
CPAC=optimvar('CPAC','LowerBound',5,'UpperBound',15);
CPEC=optimvar('CPEC','LowerBound',5,'UpperBound',15);
CPELEZ=optimvar('CPELEZ','LowerBound',5,'UpperBound',15);
CQgshp=optimvar('CQgshp','LowerBound',5,'UpperBound',15);
CQgb=optimvar('CQgb','LowerBound',5,'UpperBound',15);

Pgrid=optimvar('Pgrid',nHours,1,'LowerBound',0,'UpperBound',5); % for grid power
TSA=optimvar('TSA',nHours,1,'LowerBound',0,'UpperBound',300); % for grid power
Psell=optimvar('Psell',nHours,1,'LowerBound',0,'UpperBound',inf); % for selling power to grid
Et=optimvar('Et',nHours,1,'LowerBound',0,'UpperBound',10);
Eq=optimvar('Eq',nHours,1,'LowerBound',0,'UpperBound',10);
Ec=optimvar('Ec',nHours,1,'LowerBound',0,'UpperBound',10);
Eh2=optimvar('Eh2',nHours,1,'LowerBound',0,'UpperBound',10);
Pbch=optimvar('PbCh',nHours,1,'LowerBound',0,'UpperBound',1); % for battery charge
Pbdch=optimvar('Pbdch',nHours,1,'LowerBound',0,'UpperBound',1); % for battery discharge
qch=optimvar('qch',nHours,1,'LowerBound',0,'UpperBound',1); % for TES charge
qdch=optimvar('qdch',nHours,1,'LowerBound',0,'UpperBound',1); % for TES discharge
h2ch=optimvar('h2Ch',nHours,1,'LowerBound',0,'UpperBound',1); % for H2S charge
h2dch=optimvar('h2dch',nHours,1,'LowerBound',0,'UpperBound',1); % for H2S discharge
cch=optimvar('cch',nHours,1,'LowerBound',0,'UpperBound',1); % for cold  storage charge
cdch=optimvar('cdch',nHours,1,'LowerBound',0,'UpperBound',1); % for cold storage discharge
NOMb=optimvar('NOMb','LowerBound',5,'UpperBound',50); % nominal capacity of storage
NOMq=optimvar('NOMq','LowerBound',5,'UpperBound',50); % nominal capacity of storage
NOMc=optimvar('NOMc','LowerBound',5,'UpperBound',50); % nominal capacity of storage
NOMh=optimvar('NOMh','LowerBound',5,'UpperBound',50); % nominal capacity of storage


% ===================Initaiate the binary variables for selection==========
cp1=optimvar('cp1','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp2=optimvar('cp2','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp3=optimvar('cp3','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp4=optimvar('cp4','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp5=optimvar('cp5','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp6=optimvar('cp6','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp7=optimvar('cp7','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of power eq
cp8=optimvar('cp8','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cp9=optimvar('cp9','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq1=optimvar('cq1','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq2=optimvar('cq2','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq3=optimvar('cq3','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
cq4=optimvar('cq4','Type','integer','LowerBound',0,'UpperBound',1); % optimal configuration of heating eq
%========================operation binary variables=============
isONbch=optimvar('isONbch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge
isONbdch=optimvar('isONbdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge
isONqch=optimvar('isONqch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); 
isONqdch=optimvar('isONqdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
isONcch=optimvar('isONcch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); 
isONcdch=optimvar('isONcdch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
isONh2ch=optimvar('isONh2Ch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); 
isONh2dch=optimvar('isONh2dch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
w1=optimvar('w1',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1); %binary to buy power
w2=optimvar('w2',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1); %binary to sell power
b1=optimvar('b1',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
b2=optimvar('b2',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
b3=optimvar('b3',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
b4=optimvar('b4',nHours,1,'Type','integer','LowerBound',0,'Upperbound',1);
npv=optimvar('npv',nHours,1,'Type','integer','LowerBound',0,'UpperBound',5);  %number of pv to be installed
nwt=optimvar('nwt',nHours,1,'Type','integer','LowerBound',0,'UpperBound',5);    %number of wt to be installed

                %==============battery model============
mesprob.Constraints.Batstoredele=selfdchb*Et(idxHr2Toend-1,:)+(echb.*Pbch(idxHr2Toend,:)*dt)-...
    ((Pbdch(idxHr2Toend,:)*dt)./edchb)==Et(idxHr2Toend,:);
mesprob.Constraints.SOCb=Et(Time(1),:)<=Et(idxHr2Toend(end),:);
mesprob.Constraints.Etout=Et<= NOMb;
                %===============thermal storage model===========
mesprob.Constraints.thermalstored=(1-qselfdch)*Eq(idxHr2Toend-1,:)+(echq.*(qch(idxHr2Toend,:)*dt))-...
    ((qdch(idxHr2Toend,:)*dt)./edchq)==Eq(idxHr2Toend,:);
mesprob.Constraints.SOCq=Eq(Time(1),:)<=Eq(idxHr2Toend(end),:);
mesprob.Constraints.Eqout=Eq<= NOMq;
                %==========cold storage model==========
mesprob.Constraints.coldstored=(1-cselfdch)*Ec(idxHr2Toend-1,:)+(echc.*(cch(idxHr2Toend,:)*dt))...
    -((cdch(idxHr2Toend,:)*dt)./edchc)==Ec(idxHr2Toend,:);
mesprob.Constraints.SOCc=Ec(Time(1),:)<=Ec(idxHr2Toend(end),:);
mesprob.Constraints.Ecout=Ec<= NOMc;

%================hydrogen storage model====================================
mesprob.Constraints.hydeogenstored=selfdchh*Eh2(idxHr2Toend-1,:)+(echh2*(h2ch(idxHr2Toend,:)*dt))...
    -((h2dch(idxHr2Toend,:)*dt)./edchh2)==Eh2(idxHr2Toend,:);
mesprob.Constraints.SOCh=Eh2(Time(1),:)<=Eh2(idxHr2Toend(end),:);
mesprob.Constraints.Ehout=Eh2 <= NOMh;

%===================all equality constraints===============================
PVt=npv*PR*fR.*(SolarIrradiation'/GS);
PWt=nwt*Pwind;                              %wind power output
Pts=TSeff*TSA.*SolarIrradiation';           %thermal collector
Prce=gLHV*grce*RCEele;                      %electric output of RCE
Pice=(gLHV*iceelec*iceg);                   %electric output of ICE
Ppemfc=hLHV*h1*Pemeff;
Psofc=((hLHV*h2.*b1)+(gLHV*GgSOF.*b2))*sofeff;
Pmofc=((hLHV*h3.*b3)+(gLHV*Gpmofc.*b4))*Mofceff;
Helz=(elez/hLHV).*Pelez;                         %electrolyser output
AC=copAC*QAC;                                    %absorption chiller output
EC=copEC*Pecin;                                %electric chiller output
Qashp=Ashpeff*Pashp;                            %air source hp
Qgshp=gbheff*gbGsh*gLHV;                        %gas source hp
Qwshp=gshpef*Pgshp;                             %ground source hp
Qgb=gbeff*gbGas*gLHV;                           %gas boiler output

%=================capacity sizing==========================================
mesprob.Constraints.consa2=PVt<=CPV; 
mesprob.Constraints.consa1=CPVmin*cp1<=CPV; mesprob.Constraints.consa=CPV<=CPVmax*cp1; 
mesprob.Constraints.consi2=PWt<=CPW; 
mesprob.Constraints.consi1=CPWmin*cp2<=CPW; mesprob.Constraints.consi=CPW<=CPWmax*cp2; 
mesprob.Constraints.consb2=Pts<=CPTS;
mesprob.Constraints.consb1=CPTSmin*cp3<=CPTS; mesprob.Constraints.consb=CPTS<=CPTSmax*cp3;
mesprob.Constraints.consj2=Prce<=CPRCE; 
mesprob.Constraints.consj1=CPRCEmin*cp4<=CPRCE; mesprob.Constraints.consj=CPRCE<=CPRCEmax*cp4; 
mesprob.Constraints.consc2=Pice<=CPICE;
mesprob.Constraints.consc1=CPICEmin*cp5<=CPICE; mesprob.Constraints.consc=CPICE<=CPICEmax*cp5;
mesprob.Constraints.consk2=Ppemfc<=CPEM; 
mesprob.Constraints.consk1=CPEMmin*cp6<=CPEM; mesprob.Constraints.consk=CPEM<=CPEMmax*cp6; 
mesprob.Constraints.consd2=Psofc<=CPSOF;
mesprob.Constraints.consd1=CPSOFmin*cp7<=CPSOF; mesprob.Constraints.consd=CPSOF<=CPSOFmax*cp7;
mesprob.Constraints.consl=Pmofc<=CPMOF; 
mesprob.Constraints.consl1=CPMOFmin*cp8<=CPMOF; mesprob.Constraints.consl=CPMOF<=CPMOFmax*cp8; 
mesprob.Constraints.conse2=Helz<=CPELEZ;
mesprob.Constraints.conse1=CPELEZmin*cp9<=CPELEZ; mesprob.Constraints.conse=CPELEZ<=CPELEZmax*cp9;
mesprob.Constraints.consn2=Qashp<=CQashp; 
mesprob.Constraints.consn1=CQashp<=CQashpmin*cq1; mesprob.Constraints.consn=CQashp<=CQashpmax*cq1; 
mesprob.Constraints.consg2=Qgshp<=CQgshp;
mesprob.Constraints.consg1=CQgshpmin*cq2<=CQgshp; mesprob.Constraints.consg=CQgshp<=CQgshpmax*cq2;
mesprob.Constraints.consh2=Qgb<=CQgb;
mesprob.Constraints.consh1=CQgbmin*cq3<=CQgb; mesprob.Constraints.consh=CQgb<=CQgbmax*cq3;
mesprob.Constraints.conso2=Qwshp<=CQwshp; 
mesprob.Constraints.conso1=CQwshpmin*cq4<=CQwshp; mesprob.Constraints.conso=CQwshp<=CQwshpmax*cq4; 
mesprob.Constraints.consm=AC<=CPAC; 
mesprob.Constraints.consf=EC<=CPEC; 

Qpemfc=Ppemfc*((1-Pemeff-phloss)/Pemeff);
Qsofc=Psofc*((1-sofeff-sofhloss)/sofeff);
Qmofc=Pmofc*((1-Mofceff-mofhloss)/Mofceff);
Qrce=Prce*((1-RCEele-RCEloss)/RCEele);           %heat output of RCE
Qice=Pice*((1-iceelec-iceloss)/iceelec);         %heat output of ICE
Qfc=Qpemfc+Qsofc+Qmofc;                           %total heat from FC
Qchp=Qrce+Qice;                                    %total heat output from gas turbine
NG=Gpmofc+GgSOF+gbGas+gbGsh;                     %total natural gas demand

%====================all inequality constraints============================
mesprob.Constraints.econs1=Pbch<=PbChmax*isONbch; % Battery Limit of allowed Charging 
mesprob.Constraints.econs2=Pbdch<=PbDischmax*isONbdch; % Limit of allowed Discharging
mesprob.Constraints.econs3=isONbch+isONbdch<=1; % Forbidden the battery's charging and discharging simultaneously
mesprob.Constraints.econs4=qch<=qchmax*isONqch; % TES Limit of allowed Charging 
mesprob.Constraints.econs5=qdch<=qdchmax*isONqdch; % Limit of allowed Discharging
mesprob.Constraints.econs6=isONqch+isONqdch<=1;
mesprob.Constraints.ccons7=cch<=cchmax*isONcch; % CS Limit of allowed Charging 
mesprob.Constraints.econs8=cdch<=cdchmax*isONcdch; % CS Limit of allowed Discharging
mesprob.Constraints.econs9=isONcch+isONcdch<=1;
mesprob.Constraints.econs10=h2ch<=h2chmax*isONh2ch; % H2 storage Limit of allowed Charging 
mesprob.Constraints.econs11=h2dch<=h2dchmax*isONh2dch; % Limit of allowed Discharging
mesprob.Constraints.econs12=isONh2ch+isONh2dch<=1;
mesprob.Constraints.econs13=w1+w2<=1; %forbidden buy and selling same time.
mesprob.Constraints.Bsel=b1+b2<=1;
mesprob.Constraints.Bse2=b3+b4<=1;


%===============Carbon emission constraint==============
mesprob.Constraints.Cco2=sum((Gpmofc+GgSOF+gbGas+gbGsh)*dt*CO2g,1)+sum((Pgrid*dt)*CO2grid,1)<=...
    v.*((sum(LD'+HD'+CD')./Tloss).*1/grideff); %where v varies btw 0 and 1

%===============power balance=================

mesprob.Constraints.powerbalance=PVt(Time,:)+PWt(Time,:)+Prce(Time,:)+Pice(Time,:)...
    +Psofc(Time,:)+Pmofc(Time,:)+Ppemfc(Time,:)-Pbch(Time,:)+Pbdch(Time,:)+(w1.*Pgrid(Time,:))-(w2.*Psell(Time,:))==Pecin(Time,:)+Pashp(Time,:)+Pgshp(Time,:)+Pelez(Time,:)+LD';
%=============Heat balance==================
mesprob.Constraints.heatbalance= Qfc(Time,:)+Qchp(Time,:)+Qashp(Time,:)+Qwshp(Time,:)+Qgb(Time,:)+Qgshp(Time,:)+Pts(Time,:)-qch+qdch-HD'==0;
%==============hydrogen balance==============
mesprob.Constraints.hydrogenbalance=h1+h2+h3-Helz(Time,:)-h2ch+h2dch==0;
%==============cold balance==================
mesprob.Constraints.coldbalance=AC(Time,:)+EC(Time,:)-cch+cdch-CD'==0;

%=================objective function formulation================
%the total cost consist of annnualised investment cost, operation, maintenance and
%replacement cost% %annualised investment cost of each equipment
%==============investment cost of each equipment============
Cpvinv=pvcost*CPV; Cpwinv=pwcost*CPW; Cprceinv=rcecost*CPRCE; CpPTSinv=CPTS*TScost;
Cpiceinv=icecost*CPICE; Cpeminv=pemcost*CPEM; Cpsofinv=sofcost*CPSOF; Cpmofinv=mofcost*CPMOF;
Cpelzinv=elzcost*CPELEZ; CpBinv=BEcost*NOMb;
Cptesinv=TEScost*NOMq; Cpashinv=ashcost*CQashp; Cpwshpinv=wshpcost*CQwshp; Cpgbinv=gbcost*CQgb;
Cpgshpinv=gshpcost*CQgshp; CpACinv=CPAC*ACcost; CpECinv=CPEC*ECcost;
CpHSinv=HScost*NOMh; CpCSinv=NOMc*CScost;

Cinv=(Cpvinv+Cpwinv+Cprceinv+Cpiceinv+Cpeminv+Cpsofinv+Cpmofinv+Cpelzinv...
    +CpBinv+Cptesinv+Cpashinv+Cpwshpinv+Cpgbinv+Cpgshpinv+CpACinv+CpECinv+CpHSinv+CpCSinv+CpPTSinv)*AF;

%==============Maintenance and operational cost of equipment=============
Cm1=(Cprceinv*rcemc)+(Cpiceinv*icemc)+(Cpvinv*pvmc)+(Cpwinv*pwmc)+(Cpeminv*pemmc)+(Cpsofinv*sofmc)+(Cpmofinv*mofmc)...
    +(Cpelzinv+elzmc)+(CpPTSinv*TSmc)+(CpBinv*BEmc);
Cm2=(Cptesinv*TESmc)+(Cpashinv*ashmc)+(Cpwshpinv*wshpmc)+(Cpgbinv*gbmc)+(CpACinv*ACmc)+(CpECinv*ECmc)...
    +(CpHSinv*HSmc)+(CpCSinv*CSmc)+(Cpgshpinv*gshpmc);
   Cmain=Cm1+Cm2;

%===============Operation cost on fuel and gas==================
NGcost=sum(((NG*dt).*Cgas),1); gridcost=sum(((Pgrid*dt).*Cgrid'),1); sellcost=sum(((Psell*dt).*SellingCost),1);

%==============energy storages degradation cost=============
EESc1=sum(((Pbdch+Pbch)*dt)*BEdgcost,1)+sum(((qch+qdch)*dt)*TESdgcost,1);
EESc2=sum(((cch+cdch)*dt)*CSdgcost,1)+sum(((h2ch+h2dch)*dt)*HSdgcost,1);
Cop=(NGcost+gridcost+EESc1+EESc2-sellcost); %annual total operation cost

%===================replacement Cost========================
rep1=(Cprceinv*rcerep)+(Cpiceinv*icerep)+(Cpvinv*pvrep)+(Cpwinv*wtrep)+(Cpeminv*pemrep)+(Cpsofinv*sofrep)...
    +(Cpmofinv*mofrep)+(Cpelzinv*elzrep)+(CpBinv*Brep);
rep2=(Cptesinv*TESrep)+(Cpashinv*ashprep)+( Cpwshpinv*wshprep)+(Cpgbinv*gbrep)+(Cpgshpinv*gshphrep)...
    +(CpACinv*ACrep)+(CpECinv*ECrep)+(CpHSinv*HSrep)+(CpCSinv*CSrep)+(CpPTSinv*TSrep);
Creplacement=(rep1+rep2)*AF;

%===========set the objective function================
Totalcost=Cinv+Cmain+Creplacement+Cop;
mesprob.Objective=Totalcost;

%===========options for the optimization algorithm====================

% options for the optimization algorithm, here we set the max time it can run for
options = optimoptions('intlinprog','MaxTime',100000);
% call the optimization solver to find the best solution
[sol,TotalCost,exitflag,output] = solve(mesprob,'options',options);

[Matt J]

When I run the code, I'm prompted for input. I don't know what to answer.

v=input('please input the carbon emission target');
Ta=input('please input the projected service life of MES equipment');

[Walter Roberson]

Above, the poster said v=1; Ta=25.

Answer 1

[Walter Roberson]

You have nonlinear constraints but intlinprog does not accept nonlinear constraints.

The issue is with

Psofc=((hLHV*h2.*b1)+(gLHV*GgSOF.*b2))*sofeff;

where h2 and b1 are both optimization variables, so Psofc becomes a Quadratic optimization expression instead of a Linear optimization expression.

When you use intlinprog, even the constraints must be linear.