• 由于题目特殊性实现的不是太好，但是在无不等式约束的情况下可以得到理想结果，日后改进
  

  文章目录

  • • 实验二： 约束优化
    • • 实验目的及要求：
      • 实验原理：
      • • • 算法1：内点法
          • 算法2：外点法
      • 实验内容（方法和步骤） ：
      • • • 题目1 编写程序实现外点法（包含源代码和运行截图）
          • • 外点法函数源代码
            • 运行截图
          • 题目2 编写程序实现内点法（包含源代码和运行截图）
          • • 内点法函数源代码
            • 运行截图

实验二： 约束优化

实验目的及要求：

1. 熟练掌握内点法和外点法的思想和原理；

2. 编写程序实现内点法和外点法求解极值实例；

实验原理：

在$Matlab $环境下，按照要求编写函数和程序，求解实例，直至取得正确的运行结果。（写出每种算法的步骤）

算法1：内点法

1. 给定严格内点$x_0$为初始点，初始障碍因子$r_1=1$，给定误差$\epsilon >0$，并令$k=1$。
2. 构造障碍函数（选择其一）：
   • $G\left(x,r_k\right)=f(x)+r_k{\sum }_{j=1}^l\frac{1}{g_j(x)}$.
   • $G\left(x,r_k\right)=f(x)-r_k{\sum }_{j=1}^l\ln g_j(x)$.
3. 求障碍函数无约束极小$G(x,r_k)$。
4. 判断精度若：$r_{k} \sum_{j=1}^{l} \frac{1}{g_{j}(x)} \leq \varepsilon $或$\left|r_{k} \sum_{j=1}^{l} \ln g_{j}(x)\right| \leq \varepsilon$,则以$x_k$为原约束问题近似极小解，停止迭代；否则，取$x_{k+1}<x_k$，并令$k{++}$，返回步骤3，继续迭代。

算法2：外点法

1. 取第一个罚因子$M_1=0.1$ , 允许误差$\epsilon >0$，令$k=1$。

2. 求解罚函数$P(x,M)$的无约束极小问题，设极小值为$x_k$。
   $$P(x,M)=f(x)+M\left\{\sum _{i=1}^m{\left[h_i(x)\right]}^2+\sum _{j=1}^l{\left[\min \left(0,g_j(x)\right)\right]}^2\right\}$$

3. 判断精度. 若在$x_k$点，罚项$<\epsilon$，则停止计算，得到原问题 的近似极小点$x_k$ ; 反之，若存在某一个$j(1\le j\le l)$，有$-g_j\left(x_k\right)>\epsilon$或存在某$i (1\leq i\leq m ) $，有$\left|{h}{i}\left({x}{k}\right)\right|>\varepsilon$则令$M_{k+1} =CM_k$ , (例如 , 可取C= 5 , 10) , 并令$k++$，返回步骤2 。

实验内容（方法和步骤） ：

题目1 编写程序实现外点法（包含源代码和运行截图）

利用$Matlab$ 编写外点算法，实现下列优化问题的求解

$$\{\begin{array}{l}\min f(x)=x_1^3+x_2^2\\ \text{s.t. }g(x)=x_1+x_2-1\end{array}$$

外点法函数源代码

function [X1, X2] = EPM(f, h)
%       EPM     外点法
%       f       目标函数
%       h       等式约束
M = 0.1;        %M值
C = 10;         %倍率
eps =  0.001;    %误差
syms x1 x2;             %定义变量
H = matlabFunction(h*h);%函数句柄

while true
    P = f + M*(h*h);        %构造函数

    P_x1 = diff(P,x1);      %求x1偏导
    P_x2 = diff(P,x2);      %求x2偏导
    X = solve(P_x1==0,P_x2==0,x1,x2); %解方程组
    X1 = double(X.x1);
    X2 = double(X.x2);
    if M * abs(H(X1, X2)) < eps  %判断终止条件
        break;
    end
    M = M * C;              %加倍
end
%%%%%%%%%%%% 多个最优解则选择罚项最小的
i = 1;
answer = 10000000;
X = [X1 X2];
[n, ~] = size(X);
for j = 1 : 1 : n
    if abs(H(X(j, 1), X(j, 2))) < answer
        answer = abs(H(X(j, 1), X(j, 2)));
        i = j;
    end
end
%%%%%%%%%%% 更新答案
X1 = X1(i);
X2 = X2(i);
end

运行截图

• 输入

>> syms x1 x2
>> f = x1 ^ 3 + x2 ^ 2;
>> h = x1 + x2 - 1;
>> [X1,X2]=EPM(f,h)

• 输出

X1 =
    0.5486
X2 =
    0.4514

题目2 编写程序实现内点法（包含源代码和运行截图）

利用$Matlab$编写内点算法，实现下列优化问题的求解

$$\{\begin{array}{l}\min f(x)=x_1^3+x_2^2\\ \text{ s.t. }g(x)=x_1+x_2-1\end{array}$$

内点法函数源代码

function [X1, X2] = IPM(f, h)
%       IPM     内点法
%       f       目标函数
%       h       等式约束
R = 1;          %障碍因子R值
C = 0.1;        %倍率
eps =  0.001;   %误差
syms x1 x2;             %定义变量

H = matlabFunction(log(h));%函数句柄
F = matlabFunction(f);%函数句柄
while true
    G = f - R * log(h);    %构造函数
    G_x1 = diff(G,x1);      %求x1偏导
    G_x2 = diff(G,x2);      %求x2偏导
    X = solve(G_x1==0,G_x2==0,x1,x2); %解方程组
    X1 = real(double(X.x1));
    X2 = real(double(X.x2));
    if abs(R * H(X1, X2)) <= eps     % 判断终止条件
        break
    end
    R = R * C;
end
%%%%%%%%%%%% 多个最优解则选择目标函数最小的
i = 1;
answer = 10000000;
X = [X1 X2];
[n, ~] = size(X);
for j = 1:1:n
    if F(X(j,1),X(j,2)) < answer && X(j,1)>0 && X(j,2)>0
        answer = F(X(j,1),X(j,2));
        i = j;
    end
end
%%%%%%%%%%% 更新答案
X1 = X1(i);
X2 = X2(i);
end

运行截图

• 输入

>> syms x1 x2
>> f = x1 ^ 3 + x2 ^ 2;
>> h = x1 + x2 - 1;
>> [X1,X2]=IPM(f,h)

• 输出

X1 =
    0.5486
X2 =
    0.4514

• 运行截图

运行截图

• 输入

>> syms x1 x2
>> f = x1 ^ 3 + x2 ^ 2;
>> h = x1 + x2 - 1;
>> [X1,X2]=IPM(f,h)

• 输出

X1 =
    0.5486
X2 =
    0.4514