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

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=-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,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), 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 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 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=xstatus 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(x0xlb) errmsg=x0 must be in the range xlb 4 & isempty(xu) if all(size(xu)xub) errmsg=x0 must be in the range x0 5 if isempty(A) & size(A,2)=size(x0,1), errmsg=Matrix A not correct.; return; end;else A=; end;if nargin6 if 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), errmsg=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 nargin8 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(size(Aeq,1),1); end;if nargin10 & 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 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=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 error.; 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)-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,:)=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(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(x0xlb)=xlb(find(x0xub)=xub(find(x0xub); 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,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); 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 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 dzsepcdzsep dzsep=dzsepc; 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(:,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); 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(part