優化之外罰函數法

罰函數法的基本思想是借助罰函數把約束問題轉化為無約束問題，然后用無約束最優方法來求解。

構造罰函數：在可行點，輔助函數的值等於原來的目標函數值；在不可行點，輔助函數值等於原來的目標函數值加上一個很大的正數。可寫成形如下式：

目標函數：

約束條件：

 

其相關代碼如下：

clc
syms x1 x2 e;                                           % e為罰因子
m(1)=1;c=10;a(1)=0;b(1)=0;                               % c為遞增系數 賦初值
f=x1^2+x2^2+e*(1-x1)^2;                                 % 構造罰函數
f0(1)=0;

%求偏導、海森陣
fx1=diff(f,'x1');
fx2=diff(f,'x2');
fx1x1=diff(fx1,'x1');
fx1x2=diff(fx1,'x2');
fx2x1=diff(fx2,'x1');
fx2x2=diff(fx2,'x2');
for k=1:100                                              %外點法e迭代循環
    x1=a(k);x2=b(k);e=m(k);
    for n=1:100                                          %牛頓法求最優值
        f1=subs(fx1);                                    %求梯度值和海森矩陣
        f2=subs(fx2);
        f11=subs(fx1x1);
        f12=subs(fx1x2);
        f21=subs(fx2x1);
        f22=subs(fx2x2);
        if(double(sqrt(f1^2+f2^2))<=0.000001)              %最優值收斂條件
            a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));
            break;
        else
            X=[x1 x2]'-inv([f11 f12;f21 f22])*[f1 f2]';
            x1=X(1,1);x2=X(2,1);
        end
    end
    if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=0.000001)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=0.000001)   %迭代收斂條件
      disp('最優坐標 x1:'),disp(a(k+1))%輸出最優點坐標，迭代次數，最優值
      disp('最優坐標 x2:'),disp(b(k+1))  
      disp('迭代次數'),disp(k)
      disp('最優值'),disp(f0(k+1))
      break;
    else
      m(k+1)=c*m(k);
    end
end

　運行結果如下：