系统优化与调度——非线性规划问题：梯度投影法之MATLAB实现

System Optimization and Scheduling 

Midterm Project 

1.       Give a non-quadratic nonlinear objective function and linear inequality constraints (at least two).

2.       Select a computational method (primal or dual) for constrained optimization you learn in the class.

3.       Show your computational results and corresponding analysis. 

4.       Submit the typed report.

5.       Attach your program. 

这是之前选的一门研究生的课的期中作业。用梯度投影法解决了非线性规划问题。开始找了一份百度文库的代码，但那个代码有问题要改好多。现在至少自己试的几个test case和MATLAB自带的函数计算结果相同。

main.m

syms x1 x2 x3;% f=2*x1^2+2*x2^2-2*x1*x2^3-4*x1^7-6*x2;% var=[x1,x2];% valst=[-1,-1];% A=[1 1;1 5;-1 0;0 -1];% b=[2 5 0 0]';% f=x1^3+x2^2-2*x1-4*x2+6;% var=[x1,x2];% valst=[0 0];% A=[2,-1;1,1;-1,0;0,-1];% b=[1 2 0 0]';var=[x1,x2,x3];valst=[10,10,10];f=-x1*x2*x3;A=[-1,-2,-2;1,2,2];b=[0 72]';[x,mimfval]=MinRosenGradientProjectionMethod(f,A,b,valst,var)[x2,fval]=fmincon('confun',valst,A,b)

MinRosenGradientProjectionMethod.m

function [x,minf]=MinRosenGradientProjectionMethod(f,A,b,x0,var,eps)%f is the objection function;%A is the constraint matrix;%b is the right-hand-side vector of the constraints;%x0 is the initial feasible point;%var is the vector of iden[endent variable;%eps is the precision;%x  is the value of the independent variable when the objective function is minimum;%minf is the minimum value of the objective function;format long;if nargin == 5   eps=1.0e-6;endsyms l;x0=transpose(x0);n=length(var);sz=size(A);m=sz(1);% m is the number of rows of Agf=jacobian(f,var);%calculate the jacobian matrix of the objective functionbConti=1;while bConti        k=0;    s=0;    A1=A;    A2=A;    b1=b;    b2=b;    for i=1:m        dfun=A(i,:)*x0-b(i); %separate matrix A and b        if abs(dfun)<0.0000001  %find matrixs that satisfy A1 x_k=b1           k=k+1;           A1(k,:)=A(i,:);              b1(k,1)=b(i);             else%find matrixs that satisfy A2 x_k<b2           s=s+1;           A2(s,:)=A(i,:);              b2(s,1)=b(i);               end    end    if k>0       A1=A1(1:k,:);       b1=b1(1:k,:);    end    if s>0       A2=A2(1:s,:);       b2=b2(1:s,:);    end    while 1        P=eye(n,n);        if k>0           tM=transpose(A1);           P=P-tM*inv(A1*tM)*A1; %calculate P;        end        gv=Funval(gf,var,x0);        gv=transpose(gv);         d=-P*gv;   %calculate the searching direction        flg=1;        if(P==zeros(n))            flg =0;          end            if flg==1            d=d/norm(d); %normorlize the searching direction        end        if d==0           if k==0              x=x0;              bConti=0;              break;            else               w=-inv(A1*tM)*A1*gv;               if w>=0                            x=x0;                  bConti=0;                  break;               else                  [u,index]=min(w);%find the negative component in w                  sA1=size(A1);                  if sA1(1)==1                     k=0;                  else                     k=sA1(2);                       A1=[A1(1:(index-1),:);A1((index+1):sA1(1),:)]; %delete corresponding row in A1                   end                end             end          else              break;          end    end   y1=x0+l*d;%new iteration variable   tmpf=Funval(f,var,y1);   bb=b2-A2*x0;   dd=A2*d;   if dd>=0      [tmpI,lm]= ForwardBackMethod(tmpf,0,0.001);   %find the searching interval     else      lm=inf;%find lambda_max      for i=1:length(dd)          if(dd(i)>0)             if bb(i)/dd(i)<lm                lm=bb(i)/dd(i);             end          end       end   end   if lm==inf       lm=1e9;   end   [xm,minf]=GoldenSectionSearch(tmpf,0,lm,1.0e-14);  %guarantee lambda>0   %find the minimizer by one dimension searching method   tol=norm(xm*d);   if tol<eps      x=x0;      break;   end   x0=x0+xm*d;   disp('x0');   x0endminf=Funval(f,var,x);

ForwardBackMethod.m 进退法确定搜索区间

function [left,right]=ForwardBackMethod(f,x0,step)if nargin==2    step = 0.01endif nargin==1    x0=0;    step = 0.01endf0 =subs(f,findsym(f),{x0});x1=x0+step;f1=subs(f,findsym(f),{x1});if(f1<=f0)    step2=2*step;    x2=x1+step;    f2=subs(f,findsym(f),{x2});    while(f1>f2)        x0=x1;        x1=x2;        f0=f1;        f1=f2;        step=2*step;        x2=x1+step;        f2=subs(f,findsym(f),{x2});    end    left=min(x0, x2);    right=max(x0, x2);else    step=2*step;    x2=x1-step;    f2=subs(f,findsym(f),{x2});    while(f0>f2)        x1=x0;        x0=x2;        f1=f0;        f0=f2;        step=2*step;        x2=x1-step;        f2=subs(f,findsym(f),{x2});    end    left=min(x1,x2);%left end point    right=max(x1,x2);%right end pointend

GoldenSectionSearch.m 0.618搜索法确定局部最优值

function [x,minf]=GoldenSectionSearch(f,a,b,eps) %0.618 search method  to find minimum value of function fformat long;if nargin==3    eps=1.0e-6;end l=a+0.382*(b-a);u=a+0.618*(b-a);k=1;tol=b-a; while tol>eps&&k<100000    fl=subs(f ,findsym(f),l);    fu=subs(f ,findsym(f),u);    if fl>fu        a=l;        l=u;        u=a+0.618*(b-a);    else        b=u;        u=l;        l=a+0.382*(b-a);    end    k=k+1;    tol=abs(b-a);endif k==100000    disp('找不到最小值！');    x=NaN;    minf=NaN;    return;endx=(a+b)/2;%return the minimizerminf=subs(f, findsym(f),x);format short;

Funval.m

function fv=Funval(f,varvec,varval) var=findsym(f); varc=findsym(varvec); s1=length(var); s2=length(varc); m=floor((s1-1)/3+1); varv=zeros(1,m); if s1 ~= s2%if the number of variable is different, deal with it specially    for i=0: ((s1-1)/3)        k=findstr(varc,var(3*i+1));        index=(k-1)/3;        %varv(i+1) = varval(index+1);        index(i+1);        varv(i+1)=varval(index(i+1));    end    fv=subs(f,var,varv);else    fv=subs(f,varvec,varval); %return the value of the functionend% disp('here')% f% varvec% varval%fv = subs(f,varvec,varval);

confun.m 凸函数

function f=confun(x)f=-x(1)*x(2)*x(3);%f=2*x(1)^2+2*x(2)^2-2*x(1)*x(2)^3-4*x(1)^7-6*x(2);%f=x(1)^3+x(2)^2-2*x(1)-4*x(2)+6;%f=x(1)^2+x(2)^2+6*x(1)+2;end