多目标非线性规划程序(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);%MATLAB5.3Optimization toolbox 2.0?% Minimize F(x) %subject to: xlb = x =xub% A*x = B% Aeq*x=Beq% C(x)=0% x1+2*x2+2*x3=72% 10=x2=20% x1-x2=10% Mdiscfun.m% function f=discfun(x)% f=-

2、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 several possibly integer variables.% Usage: errmsg,Z,X,t,c,f

3、ail = % 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), EXITFLAG = -1; end;% In optim/private/qpsub.m ($Revision: 1.21 $ $Date: 1998/09/01 21:37:56

4、 $):% 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 exist(maxSQPiter,var), maxSQPiters = maxSQPiter; else maxSQPiters=inf; end; % I guess there was a reason to put

5、maxSQPiters at infinity, but this works fine for me.global maxSQPiter;% STEP 0 CHECKING INPUTZ=; X=; t=0; c=0; fail=0;if nargin1, errmsg=x0 must be a column vector.; return; end;xstatus=zeros(size(x0);if nargin2 & isempty(xstat)if all(size(xstat)=size(x0) xstatus(1:size(xstat)=xstat;else errmsg=xsta

6、tus must be a column vector the same size as x0.; return;end;if any(xstatus=round(xstatus) | xstatus0 | 23 & 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 4 & isem

7、pty(xu)if all(size(xu)xub) errmsg=x0 must be in the range x0 5if isempty(A) & size(A,2)=size(x0,1), errmsg=Matrix A not correct.; return; end;else A=; end;if nargin6if isempty(B) & any(size(B)=size(A,1) 1), errmsg=Column vector B not correct.; return; end;else B=; end;if isempty(A) & isempty(B), err

8、msg=A and B should only be nonempty together.; return; end;if isempty(B) & isempty(A), B=zeros(size(A,1),1); end;if nargin7 & isempty(Aeq)if size(Aeq,2)=size(x0,1), errmsg=Matrix Aeq not correct.; return; end;else Aeq=; end;if nargin8if isempty(Beq) & any(size(Beq)=size(Aeq,1) 1), errmsg=Column vect

9、or 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(size(Aeq,1),1); end;if nargin10 & isempty(setts)if all(size(setts)=size(settings) settings(setts=0)=setts(setts=0

10、);else errmsg=settings should be a row vector of length 3.; return; end;end;if nargin12, options1=; end;options1=optimset(optimset(fmincon),options1);if nargin13, options2=; end;options2=optimset(optimset(fmincon),options2);if nargin0 & round(maxSQPit)=maxSQPit maxSQPiter=maxSQPit;else errmsg=maxSQP

11、iter 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 error.; return;);if size(C,2)1 | size(Ceq,2)1, errmsg=C en Ceq must be column vectors.; return; end;end;% STEP 1 INITIA

12、LISATIONcurrentwarningstate=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)-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);st

13、ackxub(:,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,:)=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, e

14、rrmsg=x0 not correct.; return; end;end;end;end;end;end;errx=optimget(options2,TolX);errcon=optimget(options2,TolCon);fail=0;c=0;% STEP 2 TERMINIATIONwhile stacksize0 c=c+1;% STEP 3 LOADING OF CSP x0=stackx0(:,stacksize); xlb=stackxlb(:,stacksize); xub=stackxub(:,stacksize); x0(find(x0 xlb)=xlb(find(

15、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)=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

16、,xub,nonlcon,varargin:);if con0 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 conerrcon, convflag=1; end;end;else convflag=feasflag; end;% STEP 5 FATHOMING K = find(xstatus=1 & xlb=xub);

17、separation=1;if convflag=z_incumbent & convflag0% FC 2 separation=0;elseif all(abs(round(x(K)-x(K)0% FC 3 z_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=0 xsep

18、c = x; xsepc(K(i)=round(x(K(i); dzsepc = abs(feval(fun,xsepc,varargin:)-z);if dzsepcdzsep dzsep=dzsepc; ixsep=K(i);end;end;end;% STEP 7 SEPARATIONif 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(

19、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,stack

20、size)=domainB(2);if domainA(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=xch

21、oice(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(