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',
Related Question
Best Answer