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

1、function errmsg,Z,X,t,c,fail = BNB18(fun,x0,xstat,xl,xu,A,B,Aeq,Beq,nonlcon,setts,options1,options2,maxSQPit,varargin);%·ÇÏßÐÔÕûÊý¹æ»®Ä£ÐÍÇó½â·ÖÖ§¶¨½ç

2、µü´úËã·¨¡£ÔÚMATLAB5.3ÖÐʹÓã¬ÐèOptimization toolbox 2.0Ö§³Ö?% Minimize F(x) %subject to: xlb <= x <=xub% A*x <= B% Aeq*x=Beq% C(x)<=0 % Ceq(x)=0% x(i)¿&#

3、201;ΪÁ¬Ðø±äÁ¿£¬ÕûÊý£¬»ò¹Ì¶¨Öµ% ʹÓøñʽ%errmsg,Z,X=BNB18('fun',x0,xstat,xl,xu,A,B,Aeq,Beq,'nonlcon',setts)%fun

4、63;º MÎļþÃû£¬±íʾ×îС»¯Ä¿±êº¯Êýf=fun(x)%x0: ÁÐÏòÁ¿£¬±íʾ±äÁ¿³õÖµ

5、%xstat£º ÁÐÏòÁ¿£¬xstat(i)=0±íʾx(i)ΪÁ¬Ðø±äÁ¿£¬1±íʾÕûÊý£¬2±íʾ¹Ì¶¨Öµ%

6、xl£º ÁÐÏòÁ¿£¬±íʾ±äÁ¿Ï½ç%xu: ÁÐÏòÁ¿£¬±íʾ±äÁ¿ÉϽç%A: ¾ØÕó, ±

7、7;ʾÏßÐÔ²»µÈʽԼÊøϵÊý%B: ÁÐÏòÁ¿, ±íʾÏßÐÔ²»µÈʽԼÊøÉϽç%Aeq

8、: ¾ØÕó, ±íʾÏßÐÔµÈʽԼÊøϵÊý%Beg: ÁÐÏòÁ¿, ±íʾÏßÐÔ²»µÈʽԼÊø

9、ÓÒ¶ËÖµ%nonlcon: MÎļþÃû£¬±íʾ·ÇÏßÐÔÔ¼Êøº¯ÊýC,Ceq=nonlin(x),ÆäÖÐC(x)Ϊ²»µÈʽÔ

10、;¼Êø, % Ceq(x)ΪµÈʽԼÊø%setts: Ëã·¨ÉèÖÃ%errmsq: ·µ»Ø´íÎóÌáʾ%Z: ·µ»ØÄ¿±êº¯Ê

11、3;×îСֵ%X: ·µ»Ø×îÓŽâ%ÀýÌâ% max x1*x2*x3% -x1+2*x2+2*x3>=0% x1+2*x2+2*x3<=72% 10<=x2<=20% x1-x2=10% ÏÈд Mº¯Êýdiscfun.m% function f=discfun(x)% f=-

12、x(1)*x(2)*x(3);%Çó½â% 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% BNB18 Finds the constrained minimum of a function of

13、several possibly integer variables.% Usage: errmsg,Z,X,t,c,fail = % BNB18(fun,x0,xstatus,xlb,xub,A,B,Aeq,Beq,nonlcon,settings,options1,options2,maxSQPiter,P1,P2,.) % BNB solves problems of the form:% Minimize F(x) subject to: xlb <= x0 <=xub% A*x <= B Aeq*x=Beq% C(x)<=0 Ceq(x)=0% x(i) is

14、 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 minimized and should return a scalar. F(x)=fev

15、al(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 constrains.% B and Beq are column vectors for

16、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.% Z is the scalar result of the minimization, X the values of

17、 the accompanying variables. % t is the time elapsed while the algorithm BNB has run, c is the number of 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 phase 1 by relaxation. This sometimes makes

18、 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 no feasible solution BNB% will not try to find a best solution.% settings(2) (standard 0) if 1: if the sub-problem did not co

19、nverge 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 branched. If such% a problem does not converg

20、e 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 converge but did return a feasible% point will be considered convergent. This might be useful if fmincon is having a hard time with% a certain problem bu

21、t 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 1000 by default.% P1,P2,. are paramet

22、ers 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 stop fmincon from hanging make the

23、following chances:% In optim/private/nlconst.m ($Revision: 1.20 $ $Date: 1998/08/24 13:46:15 $):% Get EXITFLAG independent of verbosity.% After the lines: disp(' less than 2*options.TolFun but constraints are not satisfied.') % end% EXITFLAG = -1; % end% end% status=1;% add the line: if (str

24、ncmp(howqp, 'i',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

25、= Inf;% to the line: if exist('maxSQPiter','var'), maxSQPiters = maxSQPiter; else maxSQPiters=inf; end; % I guess there was a reason to put maxSQPiters at infinity, but this works fine for me.global maxSQPiter; % STEP 0 CHECKING INPUTZ=; X=; t=0; c=0; fail=0;if nargin<2, errmsg=&#

26、39;BNB needs at least 2 input arguments.' return; end;if isempty(fun), errmsg='No fun found.' return; end;if isempty(x0), errmsg='No x0 found.' return;elseif size(x0,2)>1, errmsg='x0 must be a column vector.' return; end;xstatus=zeros(size(x0);if nargin>2 & isem

27、pty(xstat) if all(size(xstat)<=size(x0) xstatus(1:size(xstat)=xstat; else errmsg='xstatus must be a column vector the same size as x0.' return; end; if any(xstatus=round(xstatus) | xstatus<0 | 2<xstatus) errmsg='xstatus must consist of the integers 0,1 en 2.' return; end;end

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

29、f 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)<=size(x0) xub(1:size(xu,1)

30、=xu; else errmsg='xub must be a column vector the same size as x0.' return; end;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(i) is an integer variabale.

31、' return;end;xub(find(xstatus=2)=x0(find(xstatus=2);if nargin>5 if isempty(A) & size(A,2)=size(x0,1), errmsg='Matrix A not correct.' return; end;else A=; end;if nargin>6 if isempty(B) & any(size(B)=size(A,1) 1), errmsg='Column vector B not correct.' return; end;else

32、 B=; end;if isempty(A) & isempty(B), errmsg='A and B should only be nonempty together.' return; end;if isempty(B) & isempty(A), B=zeros(size(A,1),1); end;if nargin>7 & isempty(Aeq) if size(Aeq,2)=size(x0,1), errmsg='Matrix Aeq not correct.' return; end;else Aeq=; end;i

33、f nargin>8 if isempty(Beq) & any(size(Beq)=size(Aeq,1) 1), errmsg='Column vector Beq not correct.' return; end;else Beq=; end;if isempty(Aeq) & isempty(Beq), errmsg='Aeq and Beq should only be nonempty together' return; end;if isempty(Beq) & isempty(Aeq), Beq=zeros(siz

34、e(Aeq,1),1); end;if nargin<10, nonlcon='' 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 3.' return; end;end;if nargin<12, options1=; end;op

35、tions1=optimset(optimset('fmincon'),options1);if nargin<13, options2=; end;options2=optimset(optimset('fmincon'),options2);if nargin<14, maxSQPiter=1000; elseif isnumeric(maxSQPit) & all(size(maxSQPit)=1 & maxSQPit>0 & round(maxSQPit)=maxSQPit maxSQPiter=maxSQPit

36、;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=',nonlcon,'(x0,varargin:);','errmsg=''nonlcon caused err

37、or.'' return;'); if size(C,2)>1 | size(Ceq,2)>1, errmsg='C en Ceq must be column vectors.' return; end;end; % STEP 1 INITIALISATIONcurrentwarningstate=warning;warning off;tic;lx = size(x0,1);z_incumbent=inf;x_incumbent=inf*ones(size(x0);I = ceil(sum(log2(xub(find(xstatus=1)

38、-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; for i=1:size(Aeq,1) if Beq(i)=1 & all(Aeq(i,:)=0 | Aeq(i,:)=1) J=find(Aeq(i

39、,:)=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; if sum(x0(J)=0, errmsg='x0 not correct.' return; end; end; end; end; end;end;errx=optimget(options2,'TolX');errcon=optimget(options

40、2,'TolCon');fail=0;c=0; % STEP 2 TERMINIATIONwhile stacksize>0 c=c+1; % STEP 3 LOADING OF CSP 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; % STEP 4 RELA

41、XATION % PHASE 1 con=BNBCON(x0,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin:); if abs(con)>errcon & settings(1)=0 x1 dummy feasflag=fmincon('0',x0,A,B,Aeq,Beq,xlb,xub,nonlcon,options1,varargin:); if settings(3) & feasflag=0 con=BNBCON(x1,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin:); if con<e

42、rrcon, feasflag=1; end; end; else x1=x0; feasflag=1; end; % PHASE 2 if feasflag>0 x 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,xlb,xub,nonlcon,varargin:); if con<errcon, convflag=1; end; end; else convflag

43、=feasflag; end; % STEP 5 FATHOMING K = find(xstatus=1 & xlb=xub); separation=1; if convflag<0 | (convflag=0 & settings(2) % FC 1 separation=0; elseif z>=z_incumbent & convflag>0 % FC 2 separation=0; elseif all(abs(round(x(K)-x(K)<errx) & convflag>0 % FC 3 z_incumbent =

44、 z; x_incumbent = x; separation = 0; end; % STEP 6 SELECTION if 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=0 xsepc = x; xsepc(K(i)=round(x(K(i); dzsepc = abs(feval(fun,xsepc,varargin:)-z); if dzsepc>dzsep dzsep=dzsep

45、c; ixsep=K(i); end; end; end; % STEP 7 SEPARATION if xchoice(ixsep)=0 % XCHOICE=0 branch=1; domain=xlb(ixsep) xub(ixsep); while branch=1 xboundary=(domain(1)+domain(2)/2; if x(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(:,stacks