多目標(biāo)非線性規(guī)劃程序

function[errmsg,Z,X,t,c,fail]=BNB18(fun,x0,xstat,xl,xu,A,B,Aeq,Beq,nonlcon,setts,options1,options2,maxSQPit,varargin);%-QiBB6OuEy1???A￡BIQ6^a-OO§i^^QpuzuEa-^j￡6UMATLAB5.3OBE16A￡-BeOptimizationtoolbox2.0O§3O?

%MinimizeF(x)%subjectto:xlb<=x<=xub%A*x<=B%Aeq*x=Beq%C(x)<=0%Ceq(x)=0%%x(i)iEiaA-B0±aA￡￡-OuEy￡-?diii^Op%Ei6AinE^%[errmsg,Z,X]=BNB18('fun',x0,xstat,xl,xu,A,B,Aeq,Beq,'nonlcon',setts)%fun￡°MiA^^Au￡-±iE^xiBj?—A￡±eo—Eyf=fun(x)%x0:ABidA￡￡-±iE^±aA￡3oQp%xstat￡°ABidAi￡-xstat(i)=0±iE^x(i)iaA-B0±aAi￡-1±iE^OuEy￡-2±iE^iii^Op%xl￡°ABidA￡￡-±iE^±aA￡iA^Q%xu:ABidA￡￡-±iE^±aA￡Ei^Q%A:初006,±iE^iBB62?pEE^6^E0ipEy%B:ABIdAi，±iE^iBB62?pEE^6^E0Ei^Q%Aeq:初006,±iE^iBB6pEE^6^E0ipEy%Beg:ABIdAi,±iE^iBB62?pEE^6^E06diEOp%nonlcon:

MiA^^Au￡-±iE^-QiBB66^E0°_Ey[C,Ceq]=nonlin(x),^aOBC(x)ia2?pEE^6^E0,%Ceq(x)IapEE^6^E0%setts:Ea-^EeOA%errmsq:,p?0zil6iaE^%Z:?p?0Ai±e°—EyxiBjOp%X:?p?0xi6A^a%%Ayia%maxx1*x2*x3%-x1+2*x2+2*x3>=0%x1+2*x2+2*x3<=72%10<=x2<=20%x1-x2=10%I宜B'M°—Eydiscfun.m%functionf=discfun(x)%f=-x(1)*x(2)*x(3);%Q6^a%clear;x0=[25,15,10]';xstat=[111]';%xl=[2010-10]';xu=[302020]';%A=[1-2-2;122];B=[072]';Aeq=[1-10];Beq=10;%[err,Z,X]=BNB18('discfun',x0,xstat,xl,xu,A,B,Aeq,Beq);%XMAX=X',ZMAX=-Z%%BNB18Findstheconstrainedminimumofafunctionofseveralpossiblyintegervariables.

%Usage:[errmsg,Z,X,t,c,fail]=%BNB18(fun,x0,xstatus,xlb,xub,A,B,Aeq,Beq,nonlcon,settings,options1,options2,maxSQPiter,P1,P2,...)%%BNBsolvesproblemsoftheform:

%MinimizeF(x)subjectto:xlb<=x0<=xub

%%%%%

A*x<=BAeq*x=BeqC(x)<=0Ceq(x)=0x(i)iscontinuousforxstatus(i)=0x(i)integerforxstatus(i)=1x(i)fixedforxstatus(i)=2%%BNBuses:

%OptimizationToolboxVersion2.0(R11)09-Oct-1998%Fromthistoolboxfmincon.miscalled.Formoreinfotypehelpfmincon.

%%funisthefunctiontobeminimizedandshouldreturnascalar.F(x)=feval(fun,x).

%x0isthestartingpointforx.x0shouldbeacolumnvector.

%xstatusisacolumnvectordescribingthestatusofeveryvariablex(i).

%xlbandxubarecolumnvectorswithlowerandupperboundsforx.

%AandAeqarematricesforthelinearconstrains.

%BandBeqarecolumnvectorsforthelinearconstrains.

%nonlconisthefunctionforthenonlinearconstrains.

%[C(x);Ceq(x)]=feval(nonlcon,x).BothC(x)andCeq(x)shouldbecolumnvectors.

%%errmsgisastringcontaininganerrormessageifBNBfoundanerrorintheinput.

%Zisthescalarresultoftheminimization,Xthevaluesoftheaccompanyingvariables.

%tisthetimeelapsedwhilethealgorithmBNBhasrun,cisthenumberofBNBcyclesand%failisthenumberofunsolvedleafsub-problems.

%%settingsisarowvectorwithsettingsforBNB:

%settings(1)(standard0)if1:usephase1byrelaxation.Thissometimesmakesthealgorithm%faster,becausephase1meansthealgorithmfirstchecksifthereisafeasiblesolution%forasub-problembeforetryingtofindabestsolution.IfthereisnofeasiblesolutionBNB%willnottrytofindabestsolution.

%settings(2)(standard0)if1:ifthesub-problemdidnotconvergedonotbranch.Ifasub-%problemdidnotconvergethismeansBNBdidnotfindasolutionforit.NormallyBNBwill%branchtheproblemsoitcantryagaintofindasolution.

%Asub-problemthatisaleafofthebranch-and-bound-threecannotbebranched.Ifsuch%aproblemdoesnotconvergeitwillbeconsideredunfeasibleandtheparameterfailwillbe%raisedbyone.

%settings(3)(standard0)if1:if1asub-problemthatdidnotconvergebutdidreturnafeasible%pointwillbeconsideredconvergent.Thismightbeusefuliffminconishavingahardtimewith%acertainproblembutyoudowantsomeresults.

%options1andoptions2areoptionsstructuresforphase1andphase2.

%Fordetailsabouttheoptionsstructuretypehelpoptimset.

%maxSQPiterisaglobalvariableusedbyfmincon(ifmodifiedasdescribedinbnb18.m).

%maxSQPiteris1000bydefault.

%P1,P2,...areparameterstobepassedtofunandnonlcon.

%F(x)=feval(fun,x,P1,P2,...).[C(x);Ceq(x)]=feval(nonlcon,x,P1,P2,...).

%TypeeditBNB18formoreinfo.

%E.C.Kuipers%e-mail

E.C.Kuipers@cpedu.rug.nl

%FI-Lab%AppliedPhysics%RijksuniversiteitGroningen%Togetridofbugsandtostopfminconfromhangingmakethefollowingchances:%%Inoptim/private/nlconst.m($Revision:1.20$$Date:1998/08/2413:46:15$):

%GetEXITFLAGindependentofverbosity.

%Afterthelines:disp('lessthan2*options.TolFunbutconstraintsarenotsatisfied.')%end%EXITFLAG=-1;%end%end%status=1;%addtheline:if(strncmp(howqp,'i',1)&mg>0),EXITFLAG=-1;end;%%Inoptim/private/qpsub.m($Revision:1.21$$Date:1998/09/0121:37:56$):

%Stopqpsubfromhanging.

%Aftertheline:%AndyGrace7-9-90.MaryAnnBranch9-30-96.

%addtheline:globalmaxSQPiter;%andchangedtheline:maxSQPiters=Inf;%totheline:ifexist('maxSQPiter','var'),maxSQPiters=maxSQPiter;elsemaxSQPiters=inf;end;%IguesstherewasareasontoputmaxSQPitersatinfinity,butthisworksfineforme.

globalmaxSQPiter;%STEP0CHECKINGINPUTZ=[];X=[];t=0;c=0;fail=0;ifnargin<2,errmsg='BNBneedsatleast2inputarguments.';return;end;ifisempty(fun),errmsg='Nofunfound.';return;end;ifisempty(x0),errmsg='Nox0found.';return;elseifsize(x0,2)>1,errmsg='x0mustbeacolumnvector.';return;end;xstatus=zeros(size(x0));ifnargin>2&~isempty(xstat)

ifall(size(xstat)<=size(x0))

xstatus(1:size(xstat))=xstat;elseerrmsg='xstatusmustbeacolumnvectorthesamesizeasx0.';return;end;ifany(xstatus~=round(xstatus)|xstatus<0|2<xstatus)errmsg='xstatusmustconsistoftheintegers0,1en2.';return;

end;end;xlb=zeros(size(x0));xlb(find(xstatus==0))=-inf;ifnargin>3&~isempty(xl)

ifall(size(xl)<=size(x0))

xlb(1:size(xl,1))=xl;

elseerrmsg='xlbmustbeacolumnvectorthesamesizeasx0.';return;end;end;ifany(x0<xlb)

errmsg='x0mustbeintherangexlb<=x0.';return;elseifany(xstatus==1&(~isfinite(xlb)|xlb~=round(xlb)))

errmsg='xlb(i)mustbeanintegerifx(i)isanintegervariabele.';return;end;xlb(find(xstatus==2))=x0(find(xstatus==2));xub=ones(size(x0));xub(find(xstatus==0))=inf;ifnargin>4&~isempty(xu)

ifall(size(xu)<=size(x0))

xub(1:size(xu,1))=xu;

elseerrmsg='xubmustbeacolumnvectorthesamesizeasx0.';return;end;end;ifany(x0>xub)

errmsg='x0mustbeintherangex0<=xub.';return;elseifany(xstatus==1&(~isfinite(xub)|xub~=round(xub)))

errmsg='xub(i)mustbeanintegerifx(i)isanintegervariabale.';return;end;xub(find(xstatus==2))=x0(find(xstatus==2));ifnargin>5

if~isempty(A)&size(A,2)~=size(x0,1),errmsg='MatrixAnotcorrect.';return;end;elseA=[];end;ifnargin>6

if~isempty(B)&any(size(B)~=[size(A,1)1]),errmsg='ColumnvectorBnotcorrect.';return;end;elseB=[];end;ifisempty(A)&~isempty(B),errmsg='AandBshouldonlybenonemptytogether.';return;end;ifisempty(B)&~isempty(A),B=zeros(size(A,1),1);end;ifnargin>7&~isempty(Aeq)

ifsize(Aeq,2)~=size(x0,1),errmsg='MatrixAeqnotcorrect.';return;end;elseAeq=[];end;ifnargin>8

if~isempty(Beq)&any(size(Beq)~=[size(Aeq,1)1]),errmsg='ColumnvectorBeqnotcorrect.';return;end;elseBeq=[];end;ifisempty(Aeq)&~isempty(Beq),errmsg='AeqandBeqshouldonlybenonemptytogether';return;end;ifisempty(Beq)&~isempty(Aeq),Beq=zeros(size(Aeq,1),1);end;ifnargin<10,nonlcon='';end;settings=[000];ifnargin>10&~isempty(setts)

ifall(size(setts)<=size(settings))

settings(setts~=0)=setts(setts~=0);

elseerrmsg='settingsshouldbearowvectoroflength3.';return;end;end;ifnargin<12,options1=[];end;options1=optimset(optimset('fmincon'),options1);ifnargin<13,options2=[];end;options2=optimset(optimset('fmincon'),options2);ifnargin<14,maxSQPiter=1000;elseifisnumeric(maxSQPit)&all(size(maxSQPit))==1&maxSQPit>0&round(maxSQPit)==maxSQPit

maxSQPiter=maxSQPit;elseerrmsg='maxSQPitermustbeaninteger>0';return;end;eval(['z=',fun,'(x0,varargin{:});'],'errmsg=''funcausederror.'';return;');if~isempty(nonlcon)

eval(['[C,Ceq]=',nonlcon,'(x0,varargin{:});'],'errmsg=''nonlconcausederror.'';return;');

ifsize(C,2)>1|size(Ceq,2)>1,errmsg='CenCeqmustbecolumnvectors.';return;end;end;%STEP1INITIALISATIONcurrentwarningstate=warning;warningoff;tic;lx=size(x0,1);z_incumbent=inf;x_incumbent=inf*ones(size(x0));I=ceil(sum(log2(xub(find(xstatus==1))-xlb(find(xstatus==1))+1))+size(find(xstatus==1),1)+1);stackx0=zeros(lx,I);stackx0(:,1)=x0;stackxlb=zeros(lx,I);stackxlb(:,1)=xlb;stackxub=zeros(lx,I);stackxub(:,1)=xub;stacksize=1;xchoice=zeros(size(x0));if~isempty(Aeq)

j=0;

fori=1:size(Aeq,1)

ifBeq(i)==1&all(Aeq(i,:)==0|Aeq(i,:)==1)J=find(Aeq(i,:)==1);ifall(xstatus(J)~=0&xchoice(J)==0&xlb(J)==0&xub(J)==1)ifall(xstatus(J)~=2)|all(x0(J(find(xstatus(J)==2)))==0)j=j+1;xchoice(J)=j;

ifsum(x0(J))==0,errmsg='x0notcorrect.';return;end;end;end;

end;

end;end;errx=optimget(options2,'TolX');errcon=optimget(options2,'TolCon');fail=0;c=0;%STEP2TERMINIATIONwhilestacksize>0

c=c+1;

%STEP3LOADINGOFCSP

x0=stackx0(:,stacksize);

xlb=stackxlb(:,stacksize);

xub=stackxub(:,stacksize);

x0(find(x0<xlb))=xlb(find(x0<xlb));x0(find(x0>xub))=xub(find(x0>xub));stacksize=stacksize-1;

%STEP4RELAXATION

%PHASE1

con=BNBCON(x0,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin{:});

ifabs(con)>errcon&settings(1)~=0

[x1dummyfeasflag]=fmincon('0',x0,A,B,Aeq,Beq,xlb,xub,nonlcon,options1,varargin{:});ifsettings(3)&feasflag==0con=BNBCON(x1,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin{:});ifcon<errcon,feasflag=1;end;

end;

elsex1=x0;feasflag=1;end;

%PHASE2

iffeasflag>0

[xzconvflag]=fmincon(fun,x1,A,B,Aeq,Beq,xlb,xub,nonlcon,options2,varargin{:});ifsettings(3)&convflag==0con=BNBCON(x,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin{:});ifcon<errcon,convflag=1;end;

end;

elseconvflag=feasflag;end;

%STEP5FATHOMING

K=find(xstatus==1&xlb~=xub);

separation=1;

ifconvflag<0|(convflag==0&settings(2))

%FC1

separation=0;

elseifz>=z_incumbent&convflag>0

%FC2

separation=0;elseifall(abs(round(x(K))-x(K))<errx)&convflag>0%FC3z_incumbent=z;x_incumbent=x;separation=0;

end;

%STEP6SELECTIONifseparation==1&~isempty(K)dzsep=-1;

fori=1:size(K,1)dxsepc=abs(round(x(K(i)))-x(K(i)));ifdxsepc>=errx|convflag==0xsepc=x;xsepc(K(i))=round(x(K(i)));dzsepc=abs(feval(fun,xsepc,varargin{:})-z);ifdzsepc>dzsepdzsep=dzsepc;ixsep=K(i);end;end;

end;

%STEP7SEPARATION

ifxchoice(ixsep)==0%XCHOICE==0branch=1;domain=[xlb(ixsep)xub(ixsep)];whilebranch==1xboundary=(domain(1)+domain(2))/2;ifx(ixsep)<xboundarydomainA=[domain(1)floor(xboundary)];domainB=[floor(xboundary+1)domain(2)];elsedomainA=[floor(xboundary+1)domain(2)];domainB=[domain(1)floor(xboundary)];end;stacksize=stacksize+1;stackx0(:,stacksize)=x;stackxlb(:,stacksize)=xlb;stackxlb(ixsep,stacksize)=domainB(1);stackxub(:,stacksize)=xub;stackxub(ixsep,stacksize)=domainB(2);ifdomainA(1)==domainA(2)stacksize=stacksize+1;stackx0(:,stacksize)=x;stackxlb(:,stacksize)=xlb;stackxlb(ixsep,stacksize)=domainA(1);stackxub(:,stacksize)=xub;stackxub(ixsep,stacksize)=domainA(2);branch=0;elsedomain=domainA;branch=1;end;end;

else%XCHOICE~=0L=find(