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

function[errmsg,Z,X,t,c,fail]=BNB18(fun,x0,xstat,xl,xu,A,B,Aeq,Beq,nonlcon,setts,options1,options2,maxSQPit,varargin);

%·???D???êy1????￡Dí?ó?a·??§?¨??μü′ú??·¨?￡?úMATLAB5.3?Dê1ó?￡?DèOptimizationtoolbox2.0?§3??

%MinimizeF(x)

%subjectto:xlb<=x<=xub

%A*x<=B

%Aeq*x=Beq

%C(x)<=0

%Ceq(x)=0

%

%x(i)?é?aá?D?±?á?￡???êy￡??ò1ì?¨?μ

%ê1ó???ê?

%[errmsg,Z,X]=BNB18('fun',x0,xstat,xl,xu,A,B,Aeq,Beq,'nonlcon',setts)

%fun￡oM???t??￡?±íê?×?D??ˉ??±êoˉêyf=fun(x)

%x0:áD?òá?￡?±íê?±?á?3??μ

%xstat￡oáD?òá?￡?xstat(i)=0±íê?x(i)?aá?D?±?á?￡?1±íê???êy￡?2±íê?1ì?¨?μ

%xl￡oáD?òá?￡?±íê?±?á?????

%xu:áD?òá?￡?±íê?±?á?é???

%A:???ó,±íê???D?2?μèê???ê??μêy

%B:áD?òá?,±íê???D?2?μèê???ê?é???

%Aeq:???ó,±íê???D?μèê???ê??μêy

%Beg:áD?òá?,±íê???D?2?μèê???ê?óò???μ

%nonlcon:M???t??￡?±íê?·???D???ê?oˉêy[C,Ceq]=nonlin(x),???DC(x)?a2?μèê???ê?,

%Ceq(x)?aμèê???ê?

%setts:??·¨éè??

%errmsq:·μ??′í?óìáê?

%Z:·μ????±êoˉêy×?D??μ

%X:·μ??×?ó??a

%

%àyìa

%maxx1*x2*x3

%-x1+2*x2+2*x3>=0

%x1+2*x2+2*x3<=72

%10<=x2<=20

%x1-x2=10

%?èD′Moˉêydiscfun.m

%functionf=discfun(x)

%f=-x(1)*x(2)*x(3);

%?ó?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=Beq

%C(x)<=0Ceq(x)=0

%x(i)iscontinuousforxstatus(i)=0

%x(i)integerforxstatus(i)=1

%x(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

%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;

%STEP0CHECKINGINPUT

Z=[];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;

%STEP1INITIALISATION

currentwarningstate=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;

%STEP2TERMINIATION

whilestacksize>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==0

con=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==0

con=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

%FC3

z_incumbent=z;

x_incumbent=x;

separation=0;

end;

%STEP6SELECTION

ifseparation==1&~isempty(K)

dzsep=-1;

fori=1:size(K,1)

dxsepc=abs(round(x(K(i)))-x(K(i)));

ifdxsepc>=errx|convflag==0

xsepc=x;xsepc(K(i))=round(x(K(i)));

dzsepc=abs(feval(fun,xsepc,varargin{:})-z);

ifdzsepc>dzsep

dzsep=dzsepc;

ixsep=K(i);

end;

end;

end;

%STEP7SEPARATION

ifxchoice(ixsep)==0

%XCHOICE==0

branch=1;

domain=[xlb(ixsep)xub(ixsep)];

whilebranch==1

xboundary=(domain(1)+domain(2))/2;

ifx(ixsep)<xboundary

domainA=[domain(1)floor(xboundary)];

domainB=[floor(xboundary+1)domain(2)];

else

domainA=[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;

else

domain=domainA;

branch=1;

end;

end;

else

%XCHOIC