- 由于题目特殊性实现的不是太好,但是在无不等式约束的情况下可以得到理想结果,日后改进
实验二: 约束优化
实验目的及要求:
-
熟练掌握内点法和外点法的思想和原理;
-
编写程序实现内点法和外点法求解极值实例;
实验原理:
在$Matlab $环境下,按照要求编写函数和程序,求解实例,直至取得正确的运行结果。(写出每种算法的步骤)
算法1:内点法
- 给定严格内点 x 0 x_0 x0为初始点,初始障碍因子 r 1 = 1 r_1=1 r1=1,给定误差 ε > 0 \varepsilon>0 ε>0,并令 k = 1 k=1 k=1。
- 构造障碍函数(选择其一):
- G ( x , r k ) = f ( x ) + r k ∑ j = 1 l 1 g j ( x ) G\left(x, r_{k}\right)=f(x)+r_{k} \sum_{j=1}^{l} \frac{1}{g_{j}(x)} G(x,rk)=f(x)+rk∑j=1lgj(x)1.
- G ( x , r k ) = f ( x ) − r k ∑ j = 1 l ln g j ( x ) G\left(x, r_{k}\right)=f(x)-r_{k} \sum_{j=1}^{l} \ln g_{j}(x) G(x,rk)=f(x)−rk∑j=1llngj(x).
- 求障碍函数无约束极小 G ( x , r k ) G(x,r_k) G(x,rk)。
- 判断精度若:$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:外点法
-
取第一个罚因子 M 1 = 0.1 M_1=0.1 M1=0.1 , 允许误差 ε > 0 \varepsilon>0 ε>0,令 k = 1 k=1 k=1。
-
求解罚函数 P ( x , M ) P(x,M) P(x,M)的无约束极小问题,设极小值为 x k x_k xk。
P ( x , M ) = f ( x ) + M { ∑ i = 1 m [ h i ( x ) ] 2 + ∑ j = 1 l [ min ( 0 , g j ( x ) ) ] 2 } 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\} P(x,M)=f(x)+M{i=1∑m[hi(x)]2+j=1∑l[min(0,gj(x))]2} -
判断精度. 若在 x k x_k xk点,罚项 < ε <\varepsilon <ε,则停止计算,得到原问题 的近似极小点 x k x_k xk ; 反之,若存在某一个 j ( 1 ≤ j ≤ l ) j (1\leq j \leq l ) j(1≤j≤l),有 − g j ( x k ) > ε -g_{j}\left(x_{k}\right)>\varepsilon −gj(xk)>ε或存在某$i (1\leq i\leq m ) , 有 ,有 ,有\left|{h}{i}\left({x}{k}\right)\right|>\varepsilon 则 令 则令 则令M_{k+1} =CM_k$ , (例如 , 可取C= 5 , 10) , 并令 k + + k{++} k++,返回步骤2 。
实验内容(方法和步骤) :
题目1 编写程序实现外点法(包含源代码和运行截图)
利用 M a t l a b Matlab Matlab 编写外点算法,实现下列优化问题的求解
{ min f ( x ) = x 1 3 + x 2 2 s.t. g ( x ) = x 1 + x 2 − 1 \left\{\begin{array}{l} \min f(x)=x_{1}^{3}+x_{2}^{2} \\ \text {s.t. } g(x)=x_{1}+x_{2}-1 \end{array}\right. {minf(x)=x13+x22s.t. g(x)=x1+x2−1
外点法函数源代码
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 编写程序实现内点法(包含源代码和运行截图)
利用 M a t l a b Matlab Matlab编写内点算法,实现下列优化问题的求解
{ min f ( x ) = x 1 3 + x 2 2 s.t. g ( x ) = x 1 + x 2 − 1 \left\{\begin{array}{l} \min f(x)=x_{1}^{3}+x_{2}^{2} \\ \text { s.t. } g(x)=x_{1}+x_{2}-1 \end{array}\right. {minf(x)=x13+x22 s.t. g(x)=x1+x2−1
内点法函数源代码
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