prob=optimproblem('ObjectiveSense','maximize');
%Sets
suppliers=1; 
% manufacture1=1;
% manufacture2=1;
% manufacture3=1;
% manufacture4=1;
plant=1;
retailer=1;
customer=1;
collection_centre=1;
disassembly_centre=1;
refurbishing_centre=1;
waste_centre=1;
periods=4;
%Decision variables
Qij=optimvar('Qij',suppliers,plant,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
% Qjk1=optimvar('Qjk1',manufacture2,manufacture1,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
% Qjk2=optimvar('Qjk2',manufacture2,manufacture3,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
% Qjk3=optimvar('Qjk3',manufacture4,manufacture3,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qjk=optimvar('Qjk',retailer,plant,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qjk1=optimvar('Qjk1',waste_centre,plant,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qlm=optimvar('Qlm',retailer,customer,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qmn=optimvar('Qmn',collection_centre,customer,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qno=optimvar('Qno',collection_centre,disassembly_centre,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qmj=optimvar('Qmj',collection_centre,plant,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qop=optimvar('Qop',waste_centre,disassembly_centre,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qop1=optimvar('Qop1',refurbishing_centre,disassembly_centre,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qrc=optimvar('Qrc',refurbishing_centre,plant,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
Qj=optimvar('Qj',1,plant,'Type','integer','LowerBound',0,'UpperBound',1000000000000);
njmachine=optimvar('njmachine',1,plant,'Type','integer','LowerBound',0,'UpperBound',1000000);
njcompressor=optimvar('njcompressor',1,plant,'Type','integer','LowerBound',0,'UpperBound',1000000);
nilum=optimvar('niulm',1,plant,'Type','integer','LowerBound',0,'UpperBound',1000000);
ncond=optimvar('ncond',1,plant,'Type','integer','LowerBound',0,'UpperBound',1000000);
nilumre=optimvar('nilumre',1,retailer,'Type','integer','LowerBound',0,'UpperBound',1000000);
ncondre=optimvar('ncondre',1,retailer,'Type','integer','LowerBound',0,'UpperBound',1000000);
nilumcc=optimvar('nilumcc',1,collection_centre,'Type','integer','LowerBound',0,'UpperBound',1000000);
ncondcc=optimvar('ncondcc',1,collection_centre,'Type','integer','LowerBound',0,'UpperBound',1000000);

alpha=optimvar('alpha',1,1,'Type','continuous','LowerBound',0,'UpperBound',1);
beta=optimvar('beta',1,1,'Type','continuous','LowerBound',0,'UpperBound',1);
gamma=optimvar('gamma',1,1,'Type','continuous','LowerBound',0,'UpperBound',1);
delta=optimvar('delta',1,1,'Type','continuous','LowerBound',0,'UpperBound',1);
%Objective function
prob.Objective=alpha-beta+gamma+delta;
%Constraints
%parameter
%Parameter of 2nd objective function
N_machine1=180; %installed power in kw
R1=1840; %manufacturing rate for machine kg/h
u=.6; %efficiency for machine in percentage
Energy_machine=(12000*(njmachine)*N_machine1)./(R1*u);
N_compressor=150; %power installed for compressor kw
ro_compressor=666; %capacity of compressed air m^3/h
v_compressor1=6;% compressed air used for machine i m^3/h
Energy_comp=((12000*N_compressor*v_compressor1*njcompressor)./(R1*u*ro_compressor));
G_illum=751569; %mass production per month for manufacturing system kg
N_illu=0.5; % power installed for illumation kw
Gpl=751569;
N_cond_cc=2; % power installed in collection center
G_cc=53164;
N_illumcc=.3;
Energy_illum=((12000*N_illu*nilum)./(G_illum))+((840*N_illumcc*nilumcc)./(G_cc));
G_cond=751569; %mass production per month for manufacturing system kg
N_cond=3; % power installed for illumation kw
Energy_cond=((12000*N_cond*ncond)./(G_cond))+((840*N_cond_cc*ncondcc)./(G_cc));



%parameter of first objective
cfixed=6000000; % fixed cost
craw=2; %per kg %raw material cost
cprod=3; % unit, product cost per unit
cin=2; % inventory cost per unit
v=20,000 % capacity of vehicle per kg
ctrans=2; %transportation cost
dist_sup=50; % per mile
dsp=50; %  mile distance from supplier to manufacture 1,2,3,4
dpr=30; %  mile distance from manufacture to retailer
dpw=5; %mile distance from plant to waste center
drc=15; %mile distance from retailer to customer
dcc=10;%mile distance from customer to collection center
dcd=20; %mile distance from collection center to disassembly center
ddw=8; %mile distance from disassembly center to waste center
ddr=17; %mile distance from disassembly center to refurbishing
drp=20; %mile distance from refurbashing to plant
dcp=15; %mile distance from colection centre to plant
N_re_air=0.5; %installed power in retailer
G_re=632650; %quantity storred in kg per month
N_re_ilum=0.1; %installed power for ilumination in retailer

Enery_re=(N_re_air*ncondre*12000)./(G_re)+(N_re_ilum*nilumre*12000)./(G_re);
%3rd objective paramter:
alpha_pl=0.04
alpha_t=0.41 % carbon emission during transportation
%fourvalued neutrosophic constraints
f1=cfixed+sum(sum(craw*(Qij)))+sum(sum(cprod*Qjk))+...
    sum(sum(cin*Qjk))+sum(sum(((ctrans*(Qij)*dsp)/v)))+...
sum(sum(((ctrans*(Qjk)*dpr)/v)))+sum(sum(((ctrans*(Qjk1)*dpw)/v)))+...
sum(sum(((ctrans*(Qlm)*drc)/v)))+sum(sum(((ctrans*(Qmn)*dcc)/v)))+...
sum(sum(((ctrans*(Qno)*dcd)/v)))+sum(sum(((ctrans*(Qop)*ddw)/v)))+...
sum(sum(((ctrans*(Qop1)*ddr)/v)))+sum(sum(((ctrans*(Qrc)*drp)/v)))+...
sum(sum(((ctrans*(Qmj)*dcp)/v)));
f2=Energy_machine+Energy_comp+Energy_illum+Energy_cond+Enery_re;
f3=alpha_pl*Energy_machine+Energy_comp*alpha_pl+Energy_cond*alpha_pl+Energy_illum*alpha_pl+Enery_re*alpha_pl+((alpha_t*Qij*dsp)./(v))+((alpha_t*Qjk*dpr)./(v))+((alpha_t*Qjk1*dpw)./(v))+((alpha_t*Qlm*drc)./(v))+((alpha_t*Qmn*dcc)./(v))+((alpha_t*Qmj*dcp)./(v))+((alpha_t*Qno*dcd)./(v))+((alpha_t*Qop*ddw)./(v))+((alpha_t*Qop1*ddr)./(v))+((alpha_t*Qrc*drp)./(v));
%capacity constraint
cap_con1=sum(sum(Qij))<=30000; %production capacity of supplier
cap_con2=sum(sum(Qjk))+sum(sum(Qjk1))<=16000; %capacity of plant
cap_con3=sum(Qlm)<=14000; % capacity of retailter
cap_con4=sum(Qno)+sum(Qmj)<=5000; % capacity of collection centre
cap_con5=sum(Qop)+sum(Qop1)<=3000; % capactiy of dissemabily centre
cap_con6=sum(Qrc)<=2000; %capacity of refurbishing centre
%Demand constraint
demand_1=sum(Qlm)>=12000;

rt1=0.07;% return rate from customers 
%% collection center supply constraint
Collection_center=sum(Qmn)>=rt1*sum(Qlm);
%transhipment constraint at collection center
rt2=0.02; % return rate fromm collevction center to dissasbmly 
rt3=0.05; % return rate from collection center to plant 
tarnshipment_collection1=sum(Qmn)>=rt2*sum(Qno);
tarnshipment_collection2=sum(Qmn)>=rt3*sum(Qmj);
%transhipment for disassembly
rt4=0.10; % return rate from disassembly to refurbishing
rt5=0.20; % return rerae from disassembly to disposal 
disassembly_con1=sum(Qno)>=rt4*sum(Qop1);
disassembly_con2=sum(Qno)>=rt5*sum(Qop);

balance_con1=sum(Qij)==sum(Qjk)+sum(Qjk1);
balance_con2=sum(Qjk)==sum(Qlm);
balance_con5=sum(Qlm)==sum(Qmn)
balance_con3=sum(Qmn)==sum(Qno)+sum(Qmj);
balance_con4=sum(Qno)==sum(Qop)+sum(Qop1);
%balance_con5=sum(Qop1)==sum(Qrc);
%Informational constraints
infor1=(1-.001)*Qj>=12000;
infor2=1.5*ncond>=njmachine;
infor3=nilum>=15*njmachine;
infor4=R1*njmachine>=Qj;
infor5=3*njcompressor>=njmachine;
infor6=(92)*ncondre>=1.5*Qjk;
infor7=92*nilumre>=10*Qjk;
infor8=(98)*ncondcc>=1.5*Qmn;
infor9=98*nilumcc>=10*Qmn;

% infor7=nilum>=15*njmachine;
% infor8=R1*njmachine>=Qj;
% infor9=3*njcompressor>=njmachine;


%standard constraint
prob.Constraints.cap_con1=cap_con1;
prob.Constraints.cap_con2=cap_con2
prob.Constraints.cap_con3=cap_con3;
prob.Constraints.cap_con4=cap_con4;
prob.Constraints.cap_con5=cap_con5;
prob.Constraints.cap_con6=cap_con6;
prob.Constraints.demand_1=demand_1;
prob.Constraints.Collection_center=Collection_center;
prob.Constraints.tarnshipment_collection1=tarnshipment_collection1;
prob.Constraints.tarnshipment_collection2=tarnshipment_collection2;
prob.Constraints.disassembly_con1=disassembly_con1;
prob.Constraints.disassembly_con2=disassembly_con2;
prob.Constraints.balance_con1=balance_con1;
prob.Constraints.balance_con2=balance_con2;
prob.Constraints.balance_con3=balance_con3;
prob.Constraints.balance_con4=balance_con4;
prob.Constraints.infor1=infor1;
prob.Constraints.infor2=infor2;
prob.Constraints.infor3=infor3;
prob.Constraints.infor4=infor4;
prob.Constraints.infor5=infor5;
prob.Constraints.infor6=infor6;
prob.Constraints.infor7=infor7;
prob.Constraints.infor8=infor8;
prob.Constraints.infor9=infor9;
%Membership function for f1(minimum)
Tf1=(H(1,1)-f1); % truth membership for function1
Unf1=(unU(1,1)-f1); %  uncertainty for function1
Cf1=(Uc(1,1)-f1); %  contradictory for function1
Ff1=(f1-LF(1,1));
%Membership function for f2 (minimum)
Tf2=(H(1,2)-f2)/(H(1,2)-G(1,2)); % truth membership for function3
Unf2=(unU(1,2)-f2)/(unU(1,2)-unL(1,2)); %  uncertainty for function3
Cf2=(Uc(1,2)-f2)/(Uc(1,2)-Lc(1,2)); %  contradictory for function3
Ff2=(f2-LF(1,2))/(UF(1,2)-LF(1,2)); % falsity for function3

%Membership function for f3(minimum)
Tf3=(H(1,3)-f3)/(H(1,3)-G(1,3)); % truth membership for function3
Unf3=(unU(1,3)-f3)/(unU(1,3)-unL(1,3)); %  uncertainty for function3
Cf3=(Uc(1,3)-f3)/(Uc(1,3)-Lc(1,3)); %  contradictory for function3
Ff3=(f3-LF(1,3))/(UF(1,3)-LF(1,3)); % falsity for function3

ncon1=Tf1>=0;
ncon2=Tf2>=alpha;
ncon3=Tf3>=alpha;
ncon4=Unf1>=0;
ncon5=Unf2>=gamma;
ncon6=Unf3>=gamma;
ncon7=Cf1>=0;
ncon8=Cf2>=delta;
ncon9=Cf3>=delta;
ncon10=Ff1<=0;
ncon11=Ff2<=beta;
ncon12=Ff3<=beta;
ncon13=alpha+beta+gamma+delta>=0;
ncon14=alpha>=beta;
ncon15=alpha>=gamma;
ncon16=alpha>=delta;
prob.Constraints.ncon1=ncon1;
prob.Constraints.ncon2=ncon2
prob.Constraints.ncon3=ncon3
prob.Constraints.ncon4=ncon4
prob.Constraints.ncon5=ncon5
prob.Constraints.ncon6=ncon6
prob.Constraints.ncon7=ncon7
prob.Constraints.ncon8=ncon8
prob.Constraints.ncon9=ncon9
prob.Constraints.ncon10=ncon10
prob.Constraints.ncon11=ncon11
prob.Constraints.ncon12=ncon12
prob.Constraints.ncon13=ncon13
prob.Constraints.ncon14=ncon14
prob.Constraints.ncon15=ncon15
prob.Constraints.ncon16=ncon16
[x fval]=solve(prob)
s=struct2cell(x)
Qij=cell2mat(s(1,1));
Qj=cell2mat(s(2,1));
Qjk=cell2mat(s(3,1));
Qjk1=cell2mat(s(4,1));
Qlm=cell2mat(s(5,1));
Qmj=cell2mat(s(6,1));
Qmn=cell2mat(s(7,1));
Qno=cell2mat(s(8,1));
Qop=cell2mat(s(9,1));
Qop1=cell2mat(s(10,1));
Qrc=cell2mat(s(11,1));
ncond=cell2mat(s(12,1));
ncondcc=cell2mat(s(13,1));
ncondre=cell2mat(s(14,1));
nilumcc=cell2mat(s(15,1));
nilumre=cell2mat(s(16,1));
nilum=cell2mat(s(17,1));
njcompressor=cell2mat(s(18,1));
njmachine=cell2mat(s(19,1));

f111=cfixed+sum(sum(craw*(Qij)))+sum(sum(cprod*Qjk))+...
    sum(sum(cin*Qjk))+sum(sum(((ctrans*(Qij)*dsp)/v)))+...
sum(sum(((ctrans*(Qjk)*dpr)/v)))+sum(sum(((ctrans*(Qjk1)*dpw)/v)))+...
sum(sum(((ctrans*(Qlm)*drc)/v)))+sum(sum(((ctrans*(Qmn)*dcc)/v)))+...
sum(sum(((ctrans*(Qno)*dcd)/v)))+sum(sum(((ctrans*(Qop)*ddw)/v)))+...
sum(sum(((ctrans*(Qop1)*ddr)/v)))+sum(sum(((ctrans*(Qrc)*drp)/v)))+...
sum(sum(((ctrans*(Qmj)*dcp)/v)));
f222=(12000*(njmachine)*N_machine1)./(R1*u)+((12000*N_compressor*v_compressor1*njcompressor)./(R1*u*ro_compressor))+((12000*N_illu*nilum)./(G_illum))+((840*N_illumcc*nilumcc)./(G_cc))+((12000*N_cond*ncond)./(G_cond))+((840*N_cond_cc*ncondcc)./(G_cc))+(N_re_air*ncondre*12000)./(G_re)+(N_re_ilum*nilumre*12000)./(G_re);
f333=alpha_pl*((12000*(njmachine)*N_machine1)./(R1*u)+((12000*N_compressor*v_compressor1*njcompressor)./(R1*u*ro_compressor))+((12000*N_illu*nilum)./(G_illum))+((840*N_illumcc*nilumcc)./(G_cc))+((12000*N_cond*ncond)./(G_cond))+((840*N_cond_cc*ncondcc)./(G_cc))+(N_re_air*ncondre*12000)./(G_re)+(N_re_ilum*nilumre*12000)./(G_re))+((alpha_t*Qij*dsp)./(v))+((alpha_t*Qjk*dpr)./(v))+((alpha_t*Qjk1*dpw)./(v))+((alpha_t*Qlm*drc)./(v))+((alpha_t*Qmn*dcc)./(v))+((alpha_t*Qmj*dcp)./(v))+((alpha_t*Qno*dcd)./(v))+((alpha_t*Qop*ddw)./(v))+((alpha_t*Qop1*ddr)./(v))+((alpha_t*Qrc*drp)./(v));
 disp(f111);
disp(f222);
disp(f333);