main_function.asv

这里实现了一个最优化控制的算法

%用梯度投影法求   min f(x1,x2)=(x(1)-1)^3-8*(x(1)-x(2))^2+x(1)*x(2)^2+exp(-x(1)-x(2));
%                      S.t  x(1)+2*x(2)<=12;
%                           -x1+3*x(2)<=10
%                            x1,x2>=0
%主程序
clea   
A=[-1,-2;1,-3;1,0;0,1];   
b=[-12;-10;0;0];
x=[2;5.1];    %初始点
%A=[1,-2;-1,-1;1,0;0,1];
%b=[-2;-4;0;0];
%x=[0;0];
%A=[-1,-1;1,0;0,1];
%b=[-2;0;0];
%x=[0;0];
%A=[-1,-1;-1,-5;1,0;0,1];
%b=[-2;-5;0;0];
%x=[0;0];
%A=[-1,-1;-1,-1;1,0;0,1];
%b=[-2;;-3;0;0];
%x=[0;0];

k=0;
Wucha=0.00000001;     %设定误差

mark=1;           %mark1是个标识变量，其值为0表示程序结束
times=0;          %times表示迭代次数 

while   mark==1|times<100
   kk=0;
   times=times+1
   number1=Constain_function(x,A,b) %下面是计算有效约束
   number1(1)
   [c1,c2]=size(number1);
   if number1(1)~=0
         for i=1:c2
        b1(i)=b(number1(i));
           for j=1:2
           A1(i,j)=A(number1(i),j);
       end
   end
end
   clear A;
   clear b;
   A=[-1,-2;1,-3;1,0;0,1];   
   b=[-12;-10;0;0];
   number2=Constain_function1(x,A,b,k); %下面是计算无效约束
   k=0;
   [c1,c2]=size(number2);
if number2(1)~=0
   for i=1:c2
      b2(i)=b(number2(i));
      for j=1:2
           A2(i,j)=A(number2(i),j);
      end
    end
end
   if number1(1)==0            %如果有效约束矩阵为空
         Q=eye(2);             %将Q置为单位针 
   elseif number1(1)~=0 
         M=A1;
         Q=eye(2)-M'*inv((M*M'))*M;
   end   
     
         d=-Q*Grad_function(x) %计算可行方向
     
     
     if  norm(d)<Wucha&&number1(1)==0  %如果d的模小于给定的误差，并且有效约束矩阵为空
            mark1=0;                   
            break;                     %程序停止运行，当前点就是最优点
     elseif norm(d)<Wucha&&number1(1)~=0 
            q=inv(M*M')*M*Grad_function(x);
            hh=1;                  %hh用于判断q中元素的正负
            [mm,nn]=size(q);
        for i=1:mm 
             hh=hh&&q(i)>=0;   
        end     
        if hh==1                   %hh==1表示q中所有元素都大于0
            mark1=0;               %程序结束，当前点为库恩－图克点
            break;
        else                       %下面选择q中的最小负分量
            for i=1:mm
                if q(i)<kk
                    kk=q(i);
                    k=i
                end
            end
        k=number1(k);
        [c1,c2]=size(A);           %下面去除A中对应的行
        for i=1:k-1
            bb(i)=b(i)
            for j=1:2
                AA(i,j)=A(i,j);
            end
        end
        for i=k:c1-1
            bb(i)=b(i+1)           %下面去除b中对应的行
            for j=1:2
                AA(i,j)=A(i+1,j);
            end
        end      
       
        clear A;
        clear b;
        clear A1;
        clear A2;
        clear b1;
        clear b2;
        clear M;
        A=AA;
        b=bb;
        clear AA;
        clear bb;
        continue;
 end
elseif  norm(d)>=Wucha          %当d的模大于等于给定的误差的时候执行 
    lamda_max=Lamda_function(A2,b2,d,x);   %求出最大步长
    lamda=Gold_function(x,d,lamda_max)    %求出最佳步长，使用黄金分割法
     if(lamda<0.00000001)                  %如果最佳步长过短，则程序停止执行
              mark=0;  break;
        else x=x+lamda*d;                  %计算出下一个迭代点
     end

        clear A1;
        clear A2;
        clear b1;
        clear b2;
        clear M;
        clear A;
        clear b;
A=[-1,-2;1,-3;1,0;0,1];   
b=[-12;-10;0;0];
end
end

times                 %打印出迭代次数
x                     %目标函数取极小值时的x
Aim_function(x)       %此时目标函数的值