多目标非线性规划程序(Matlab)

1、function errmsg,Z,X,t,c,fail = BNB18(fun,x0,xstat,xl,xu,A,B,Aeq,Beq,nonlcon,setts,options1,options2,ma xSQPit,varargin);% ?D?e y1?£ Di ?o ?a - ? § ?" ? uzu ? ？£ ?u MATLAB5.3?D 1 o ?£ ?De Optimizat ion toolbox 2.0?§ 3?%Minimize F(x)%subject to: xlb <= x <=xub%A*x &

2、lt;= B%Aeq*x=Beq%C(x)<=0% Ceq(x)=0%x(i)? e ?ad ?D?± ?d ?£ ?e y£ ?o 1 i ?" ?卩% e 1 o ?e?%errmsg,Z,X=BNB18('fun',x0,xstat,xl,xu,A,B,Aeq,Beq,'nonlcon',setts)%fun£o M?t? £?±ie ?x ?D?_ ?±e o_e yf=fun(x)%x0:d D?od ?£ ?±ie ?± ?d ?3?卩

3、%xstat £ o d D?od ?£ ?xstat(i)=0 ±ie ?x(i)?a d ?D?± ?d ?£ ?1 ±ie ?e y£ ?2±ie ?1 i ?" ?卩 %xl£o d D?od ?£ ?±ie ?± ?d ?%xu:d D?od ?£ ?±ie ?± ?d ?e ?%A:? o , ±ie ?D?2?ee ?e ? ey%B:d D?od ?, ±ie ?D?2?ee ?e ?e ?%Aeq:

4、 ? o , ±ie ?D?ee ?e ? ey%Beg:d D?od ?, ±ie ?D?2?ee ?e ?oo ?卩%nonlcon:M?t?£ ?±ie ? - ?D?e ?o _e yC,Ceq=nonlin(x),?DC(x)?a2?ee ?e ?,%Ceq(x)?a ee ?e?%setts: ?ee ?%errmsq:?'i ?oide?%Z:?±e o_e yx?D?%X:飞?x ?o ?a%a y ia% max x1*x2*x3% -x1+2*x2+2*x3>=0% x1+2*x2+2*x3<=72% 1

5、0<=x2<=20%x1-x2=10% ?e D' Mo_e ydiscfun.m%function f=discfun(x)%f=-x(1)*x(2)*x(3);%?o ?a%clear;x0=25,15,10'xstat=1 1 1'%xl=20 10 -10'xu=30 20 20'%A=1 -2 -2;1 2 2;B=0 72'Aeq=1 -1 0;Beq=10;% err,Z,X=BNB18('discfun',x0,xstat,xl,xu,A,B,Aeq,Beq);% XMAX=X',ZMAX=-Z

6、% BNB18 Finds the constrained minimum of a function of several possibly integer variables.% Usage: errmsg,Z,X,t,c,fail =%BNB18(fun,x0,xstatus,xlb,xub,A,B,Aeq,Beq,nonlcon,settings,options1,opti ons2,maxSQPiter,P1,P2,.)% BNB solves problems of the form:% Minimize F(x) subject to: xlb <= x0 <=xub

7、%A*x <= BAeq*x=BeqC(x)<=0 Ceq(x)=0 x(i) is continuous for xstatus(i)=0 x(i) integer for xstatus(i)= 1 x(i) fixed for xstatus(i)=2% BNB uses:% Optimization Toolbox Version 2.0 (R11) 09-Oct-1998% From this toolbox fmincon.m is called. For more info type help fmincon.% fun is the function to be m

8、inimized and should return a scalar. F(x)=feval(fun,x).% x0 is the starting point for x. x0 should be a column vector.% xstatus is a column vector describing the status of every variable x(i).% xlb and xub are column vectors with lower and upper bounds for x.% A and Aeq are matrices for the linear c

9、onstrains.% B and Beq are column vectors for the linear constrains.% nonlcon is the function for the nonlinear constrains.% C(x);Ceq(x)=feval(nonlcon,x). Both C(x) and Ceq(x) should be column vectors.% errmsg is a string containing an error message if BNB found an error in the input.X the values of

10、the accompanying% Z is the scalar result of the minimization, variables.% t is the time elapsed while the algorithm BNB has run, c is the numberof BNB cycles and% fail is the number of unsolved leaf sub-problems.% settings is a row vector with settings for BNB:% settings(1) (standard 0) if 1: use ph

11、ase 1 by relaxation. This sometimes makes the algorithm% faster, because phase 1 means the algorithm first checks if there is a feasible solution% for a sub-problem before trying to find a best solution. If there is nofeasible solution BNB% will not try to find a best solution.% settings(2) (standar

12、d 0) if 1: if the sub-problem did not converge do not branch. If a sub-% problem did not converge this means BNB did not find a solution for it.Normally BNB will% branch the problem so it can try again to find a solution.% A sub-problem that is a leaf of the branch-and-bound-three can not be branche

13、d. If such% a problem does not converge it will be considered unfeasible and the parameter fail will be% raised by one.% settings(3) (standard 0) if 1: if 1 a sub-problem that did not convergebut did return a feasible% point will be considered convergent. This might be useful if fmincon is having a

14、hard time with% a certain problem but you do want some results.% options1 and options2 are options structures for phase 1 and phase 2.% For details about the options structure type help optimset.% maxSQPiter is a global variable used by fmincon (if modified as described in bnb18.m).% maxSQPiter is 1

15、000 by default.% P1,P2,. are parameters to be passed to fun and nonlcon.% F(x)=feval(fun,x,P1,P2,.). C(x);Ceq(x)=feval(nonlcon,x,P1,P2,.).% Type edit BNB18 for more info.% E.C. Kuipers% e-mail E.C.Kuiperscpedu.rug.nl% FI-Lab% Applied Physics% Rijksuniversiteit Groningen% To get rid of bugs and to st

16、op fmincon from hanging make the following chances:% In optim/private/nlconst.m$):% Get EXITFLAG independent of verbosity.% After the lines:disp(' less than 2*options.TolFun butconstraints are not satisfied.')%end%EXITFLAG = -1;% end%end%status=1;% add the line: if (strncmp(howqp, 'i'

17、;,1) & mg > 0), EXITFLAG = -1; end;% In optim/private/qpsub.m ($Revision: 1.21 $ $Date: 1998/09/01 21:37:56$):% Stop qpsub from hanging.% After the line: % Andy Grace 7-9-90. Mary Ann Branch 9-30-96.% add the line: global maxSQPiter;% and changed the line: maxSQPiters = Inf;% to the line: if

18、exist('maxSQPiter','var'), maxSQPiters=inf; end;($Revision:1.20 $ $Date: 1998/08/24 13:46:15maxSQPitersmaxSQPiter;else% I guess there was a reason to put maxSQPiters at infinity, fine for me.global maxSQPiter;but thisworks% STEP 0 CHECKING INPUT Z=; X=; t=0; c=0; fail=0; if nargin<

19、;2, errmsg= if isempty(fun), errmsg= if isempty(x0), errmsg= elseif size(x0,2)>1, errmsg= xstatus=zeros(size(x0); if nargin>2 & isempty(xstat) if all(size(xstat)<=size(x0) xstatus(1:size(xstat)=xstat; else errmsg= 'xstatus end ;'BNB needs at least 2 input arguments.''No

20、fun found.' ; return ; end ; 'No x0 found.' ; return ;'x0 must be a column vector.'must be a column vector the same sizeif any(xstatus=round(xstatus) | xstatus<0 | 2<xstatus) errmsg= 'xstatus must consist of the integers 0,1 en 2.' end ;end ;xlb=zeros(size(x0); xlb(

21、find(xstatus=0)=-inf; if nargin>3 & isempty(xl)returnreturnas x0.' end ; end ;returnreturnif all(size(xl)<=size(x0) xlb(1:size(xl,1)=xl;else errmsg= 'xlb must be a column vector the same size as x0.' returnend ;end ;if any(x0<xlb)errmsg= 'x0 must be in the range xlb <

22、= x0.' return ;elseif any(xstatus=1 & (isfinite(xlb) | xlb=round(xlb)errmsg= 'xlb(i) must be an integer if x(i) is an integer variabele.' return ; end ;xlb(find(xstatus=2)=x0(find(xstatus=2);xub=ones(size(x0); xub(find(xstatus=0)=inf;if nargin>4 & isempty(xu)if all(size(xu)<

23、;=size(x0) xub(1:size(xu,1)=xu;else errmsg= 'xub must be a column vector the same size as x0.' returnend ;end ;if any(x0>xub)errmsg= 'x0 must be in the range x0 <=xub.' return ;elseif any(xstatus=1 & (isfinite(xub) | xub=round(xub)errmsg= 'xub(i) must be an integer if x

24、(i) is an integer variabale.' return ; end ;xub(find(xstatus=2)=x0(find(xstatus=2); if nargin>5if isempty(A) return ; end ; else A=; end ; if nargin>6if isempty(B) correct.' ; return else B=; end ;& size(A,2)=size(x0,1),& any(size(B)=size(A,1) ; end ;if isempty(A) & isempty

25、(B), errmsg= together.' ; return ; end ; if isempty(B) & isempty(A), B=zeros(size(A,1),1); if nargin>7 & isempty(Aeq)if size(Aeq,2)=size(x0,1), errmsg= end ; else Aeq=; if nargin>8end ;errmsg= 'Matrix A not1), errmsg= 'Column'A and B should only be nonemptyend ;'Mat

26、rix Aeq not correct.'if isempty(Beq) & any(size(Beq)=size(Aeq,1) 1), errmsg= vector Beq not correct.' return ; end ;else Beq=; end ; if isempty(Aeq) & isempty(Beq), together' ; return ; end ; if isempty(Beq) & isempty(Aeq), Beq=zeros(size(Aeq,1),1); if nargin<10, nonlcon=&

27、#39;' ; end ;settings = 0 0 0;if nargin>10 & isempty(setts) if all(size(setts)<=size(settings) settings(setts=0)=setts(setts=0);else errmsg= 'settings should be a row vector of length end ; if nargin<12, options1=;errmsg= 'Aeq and Beq shouldonlycorrect.'vector B notretur

28、n'Columnbe nonemptyend ;3.' ; return ; end ;end ;'fmincon' ),options1); end ;'fmincon' ),options2);options1=optimset(optimset(if nargin<13, options2=;options2=optimset(optimset(if nargin<14, maxSQPiter=1000;elseif isnumeric(maxSQPit) & all(size(maxSQPit)=1 & max

29、SQPit>0 &round(maxSQPit)=maxSQPit maxSQPiter=maxSQPit;else errmsg= 'maxSQPiter must be an integer >0' return ; end ;eval( 'z=' ,fun, '(x0,varargin:);' , 'errmsg=''fun caused error.'' return;' );if isempty(nonlcon)eval( 'C, Ceq=',nonlc

30、on,'(x0,varargin:);' , 'errmsg=''nonlconcaused error.'' return;' );vectors.'if size(C,2)>1 | size(Ceq,2)>1, errmsg= 'C en Ceq must be column return ; end ;end ;% STEP 1 INITIALISATION currentwarningstate=warning; warning off ; tic;lx = size(x0,1); z_incu

31、mbent=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

32、;for i=1:size(Aeq,1)if Beq(i)=1 & all(Aeq(i,:)=0 | Aeq(i,:)=1) J=find(Aeq(i,:)=1);if all(xstatus(J)=0 & xchoice(J)=0 & xlb(J)=0 & xub(J)=1)if all(xstatus(J)=2) | all(x0(J(find(xstatus(J)=2)=0) j=j+1; xchoice(J)=j;end ;if sum(x0(J)=0, errmsg= 'x0 not correct.' ; return end ; e

33、nd ; end ;end ;end ;errx=optimget(options2, 'TolX' ); errcon=optimget(options2, 'TolCon' );fail=0;c=0;% STEP 2 TERMINIATION while stacksize>0c=c+1;% STEP 3 LOADING OF CSP x0=stackx0(:,stacksize); xlb=stackxlb(:,stacksize); xub=stackxub(:,stacksize); x0(find(x0<xlb)=xlb(find(x0&

34、lt;xlb); x0(find(x0>xub)=xub(find(x0>xub); stacksize=stacksize-1;% STEP 4 RELAXATION% PHASE 1 con=BNBCON(x0,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin:);if abs(con)>errcon & settings(1)=0x1 dummyfeasflag=fmincon( '0' ,x0,A,B,Aeq,Beq,xlb,xub,nonlcon,options1,varargin :);if settings(3)

35、& feasflag=0 con=BNBCON(x1,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin:);if con<errcon, feasflag=1; end ;end ;else x1=x0; feasflag=1; end ;% PHASE 2if feasflag>0x z convflag=fmincon(fun,x1,A,B,Aeq,Beq,xlb,xub,nonlcon,options2,varargin :);if settings(3) & convflag=0 con=BNBCON(x,A,B,Aeq,Beq,xl

36、b,xub,nonlcon,varargin:);if con<errcon, convflag=1; end ;end ;else convflag=feasflag; end ;% STEP 5 FATHOMINGK = find(xstatus=1 & xlb=xub); separation=1;if convflag<0 | (convflag=0 & settings(2)% FC 1separation=0;elseifz>=z_incumbent & convflag>0% FC 2separation=0;elseif all(

37、abs(round(x(K)-x(K)<errx) & convflag>0% FC 3z_incumbent = z;x_incumbent = x;separation = 0;end ;% STEP 6 SELECTIONif separation = 1 & isempty(K)dzsep=-1;for i=1:size(K,1)dxsepc = abs(round(x(K(i)-x(K(i);if dxsepc>=errx | convflag=0xsepc = x; xsepc(K(i)=round(x(K(i);dzsepc = abs(feva

38、l(fun,xsepc,varargin:)-z);if dzsepc>dzsepdzsep=dzsepc; ixsep=K(i);end ;end ;end ;% STEP 7 SEPARATIONif xchoice(ixsep)=0% XCHOICE=0branch=1;domain=xlb(ixsep) xub(ixsep);while branch=1xboundary=(domain(1)+domain(2)/2;if x(ixsep)<xboundarydomainA=domain(1) floor(xboundary); domainB=floor(xboundar

39、y+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);if domainA(1)=domainA(2) stacksize=stacksize+1; stackx0(:,stacksize)=x; stackxlb(:,stacksize)=xlb; stackxlb(ixsep,stacksize)=domainA(1); stackx