（目标管理）多目标非线性
规划程序(M)
function[errmsg,Z,X,t,c,fail]=BNB18(fun,x0,xstat,xl,xu,A,B,Aeq,Beq,nonlcon,setts,opti
ons1,options2,maxSQPit,varargin);
%·ÇÏßÐÔÕûÊý¹æ»®Ä£ÐÍÇó½â·ÖÖ§¶¨½çµü´úËã·¨¡£ÔÚÖÐʹÓã¬Ðè
Ö§³Ö?
%MinimizeF(x)
%subjectto:xlb<=x<=xub
%A*x<=B
%Aeq*x=Beq
%C(x)<=0
%Ceq(x)=0
%
%x(i)¿ÉΪÁ¬Ðø±äÁ¿£¬ÕûÊý£¬»ò¹Ì¶¨Öµ
%ʹÓøñʽ
%[errmsg,Z,X]=BNB18('fun',x0,xstat,xl,xu,A,B,Aeq,Beq,'nonlcon',setts)
%fun£ºMÎļþÃû£¬±íʾ×îС»¯Ä¿±êº¯Êýf=fun(x)
%x0:ÁÐÏòÁ¿£¬±íʾ±äÁ¿³õÖµ
%xstat£ºÁÐÏòÁ¿£¬xstat(i)=0±íʾx(i)ΪÁ¬Ðø±äÁ¿£¬1±íʾÕûÊý£¬2±íʾ¹Ì¶¨Öµ
%xl£ºÁÐÏòÁ¿£¬±íʾ±äÁ¿Ï½ç
%xu:ÁÐÏòÁ¿£¬±íʾ±äÁ¿ÉϽç
%A:¾ØÕó,±íʾÏßÐÔ²»µÈʽԼÊøϵÊý
%B:ÁÐÏòÁ¿,±íʾÏßÐÔ²»µÈʽԼÊøÉϽç
%Aeq:¾ØÕó,±íʾÏßÐÔµÈʽԼÊøϵÊý
%Beg:ÁÐÏòÁ¿,±íʾÏßÐÔ²»µÈʽԼÊøÓÒ¶ËÖµ
%nonlcon:MÎļþÃû£¬±íʾ·ÇÏßÐÔÔ¼Êøº¯Êý[C,Ceq]=nonlin(x),ÆäÖÐC(x)Ϊ²»µÈʽԼÊø,
%Ceq(x)ΪµÈʽԼÊø
%setts:Ëã·¨ÉèÖÃ
%errmsq:·µ»Ø´íÎóÌáʾ
%Z:·µ»ØÄ¿±êº¯Êý×îСֵ
%X:·µ»Ø×îÓŽâ
%
%ÀýÌâ
%maxx1*x2*x3
%-x1+2*x2+2*x3>=0
%x1+2*x2+2*x3<=72
%10<=x2<=20
%x1-x2=10
%ÏÈдMº¯Êý
%functionf=discfun(x)
%f=-x(1)*x(2)*x(3);
%Çó½â
%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,P
1,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:
%(R11)09-Oct-1998
%.
%
%(x)=feval(fun,x).
%.
%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:
%faster,becausephase1meansthealgorithmfirstchecksifthereisafeasiblesolution
%
%willnottrytofindabestsolution.
%settings(2)(standard0)if1:-
%
%branchtheproblemsoitcantryagaintofindasolution.
%
%aproblemdoesnotconvergeitwillbeconsideredunfeasibleandtheparameterfailwillbe
%raisedbyone.
%settings(3)(standard0)if1:if1asub-problemthatdidnotconvergebutdidreturnafeasible
%
%acertainproblembutyoudowantsomeresults.
%options1andoptions2areoptionsstructuresforphase1andphase2.
%Fordetailsabouttheoptionsstructuretypehelpoptimset.
%maxSQPiterisaglobalvariableusedbyfmincon().
%maxSQPiteris1000bydefault.
%P1,P2,...areparameterstobepassedtofunandnonlcon.
%F(x)=feval(fun,x,P1,P2,...).[C(x);Ceq(x)]=feval(nonlcon,x,P1,P2,...).
%TypeeditBNB18formoreinfo.
%
%@
%FI-Lab
%AppliedPhysics
%RijksuniversiteitGroningen
%Togetridofbugsandtostopfminconfromhangingmakethefollowingchances:
%
%Inoptim/private/($Revision:$$Date:1998/08/2413:46:15$):
%GetEXITFLAGindependentofverbosity.
%Afterthelines:disp('lessthan2*.')
%end
%EXITFLAG=-1;
%end
%end
%status=1;
%addtheline:if(strncmp(howqp,'i',1)&mg>0),EXITFLAG=-1;end;
%
%Inoptim/private/($Revision:$$Date:1998/09/0121:37:56$):
%Stopqpsubfromhanging.
%Aftertheline:%-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.';retu
rn;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
%XCHOICE~=0
L=find(xchoice==xchoice(ixsep));
M=intersect(K,L);
[dummy,N]=sort(x(M));
part1=M(N(1:floor(size(N)/2)));part2=M(N(floor(size(N)/2)+1:size(N)));
stacksize=stacksize+1;
stackx0(:,stacksize)=x;
O=(1-sum(stackx0(part1,stacksize)))/size(part1,1);
stackx0(part1,stacksize)=stackx0(part1,stacksize)+O;
stackxlb(:,stacksize)=xlb;
stackxub(:,stacksize)=xub;
stackxub(part2,stacksize)=0;
stacksize=stacksize+1;
stackx0(:,stacksize)=x;
O=(1-sum(stackx0(part2,stacksize)))/size(part2,1);
stackx0(part2,stacksize)=stackx0(part2,stacksize)+O;
stackxlb(:,stacksize)=xlb;
stackxub(:,stacksize)=xub;
stackxub(part1,stacksize)=0;
ifsize(part2,1)==1,stackxlb(part2,stacksize)=1;end;
end;
elseifseparation==1&isempty(K)
fail=fail+1;
end;
end;
%STEP8OUTPUT
t=toc;
Z=z_incumbent;
X=x_incumbent;
errmsg='';
eval(['warning',currentwarningstate]);
functionCON=BNBCON(x,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin);
ifisempty(A),CON1=[];elseCON1=max(A*x-B,0);end;
ifisempty(Aeq),CON2=[];elseCON2=abs(Aeq*x-Beq);end;
CON3=max(xlb-x,0);
CON4=max(x-xub,0);
ifisempty(nonlcon)
CON5=[];CON6=[];
else
[CCeq]=feval(nonlcon,x,varargin{:});
CON5=max(C,0);
CON6=abs(Ceq);
end;
CON=max([CON1;CON2;CON3;CON4;CON5;CON6]);
