EDA365电子论坛网

标题: 基于matlab约束优化之惩罚函数法 [打印本页]

作者: mutougeda 时间: 2021-3-29 11:22
标题: 基于matlab约束优化之惩罚函数法

一、算法原理
1、问题引入
之前我们了解过的算法大部分都是无约束优化问题，其算法有：黄金分割法，牛顿法，拟牛顿法，共轭梯度法，单纯性法等。但在实际工程问题中，大多数优化问题都属于有约束优化问题。惩罚函数法就可以将约束优化问题转化为无约束优化问题，从而使用无约束优化算法。

2、约束优化问题的分类
约束优化问题大致分为三类：等式约束、不等式约束、等式+不等式约束。

其数学模型为：

等式约束

min f(x),x\in R^n

s.t hv(x)=0,v=1,2,…,p<n

不等式约束

min f(x),x\in R^n

s.t gu(x)\leqslant 0,u=1,2,…,p<n

等式+不等式约束问题

min f(x),x\in R^n

s.t hv(x)=0,v=1,2,…,p<n

s.t gu(x)\leqslant 0,u=1,2,…,p<n

3、惩罚函数法定义
惩罚函数法（SUMT法）又称序列无约束极小化技术，将等式约束与不等式约束的条件，经过适当定义的复合函数加到原目标函数上构造了惩罚函数，从而取消了约束，转而求解一系列无约束优化问题。

按照惩罚函数再优化过程中的迭代点是否在约束条件的可行域内，又分为内点法、外点法和混合法

内点法：迭代点再约束条件的可行域之内，只用于不等式约束。

外点法：迭代点再约束条件的可行域之外，既用于不等式约束又可用于等式约束。

4、外点惩罚函数法
等式约束：

min f=x1.2+x2.2,x\in R^n

s.t h1(x)=x1-2=0,h2(x)=x2+3=0

算法步骤

a、构造惩罚函数：F=f+M * { [ h1(x) ]^2 + [ h2(x) ]^2 } ，式中M为初始惩罚因子；

b、然后用无约束优化极值算法求解(牛顿法)；

c、如果相邻两次惩罚函数无约束最优点之间的距离足够小【norm(x1-x0)<eps】，则收敛；

否则放大惩罚因子M=C*M，式中C为罚因子放大系数；

d、转步骤a继续迭代；
/ r# l+ u2 R7 K4 H

二、源代码

%% 外点惩罚函数法-等式约束
syms x1 x2
f=x1.^2+x2.^2;
hx=[x1-2;x2+3];%列
x0=[0;0];
M=0.01;
C=2;
eps=1e-6;
[x,result]=waidian_EQ(f,x0,hx,M,C,eps)
function [x,result]=waidian_EQ(f,x0,hx,M,C,eps)
% f 目标函数
% x0 初始值
% hx 约束函数
% M 初始罚因子
% C 罚因子放大系数
% eps 容差
%计算惩罚项
CF=sum(hx.^2); %chengfa
while 1
F=matlabFunction(f+M*CF);%目标函数，使用之前的牛顿法，需要转换成句柄
x1=Min_Newton(F,x0,eps,100);
if norm(x1-x0)<eps
x=x1;
result=double(subs(f,symvar(f),x'));
break;
else
M=M*C;
x0=x1;
end
end
end
%牛顿法
function [X,result]=Min_Newton(f,x0,eps,n)
%f为目标函数
%x0为初始点
%eps为迭代精度
%n为迭代次数
TiDu=gradient(sym(f),symvar(sym(f)));% 计算出梯度表达式
Haisai=jacobian(TiDu,symvar(sym(f)));
Var_Tidu=symvar(TiDu);
Var_Haisai=symvar(Haisai);
Var_Num_Tidu=length(Var_Tidu);
Var_Num_Haisai=length(Var_Haisai);
TiDu=matlabFunction(TiDu);
flag = 0;
if Var_Num_Haisai == 0 %也就是说海塞矩阵是常数
Haisai=double((Haisai));
flag=1;
end
%求当前点梯度与海赛矩阵的逆
f_cal='f(';
TiDu_cal='TiDu(';
Haisai_cal='Haisai(';
for k=1:length(x0)
f_cal=[f_cal,'x0(',num2str(k),'),'];
for j=1: Var_Num_Tidu
if char(Var_Tidu(j)) == ['x',num2str(k)]
TiDu_cal=[TiDu_cal,'x0(',num2str(k),'),'];
end
end
for j=1:Var_Num_Haisai
if char(Var_Haisai(j)) == ['x',num2str(k)]
Haisai_cal=[Haisai_cal,'x0(',num2str(k),'),'];
end
end
end
Haisai_cal(end)=')';
TiDu_cal(end)=')';
f_cal(end)=')';
switch flag
case 0
Haisai=matlabFunction(Haisai);
dk='-eval(Haisai_cal)^(-1)*eval(TiDu_cal)';
case 1
dk='-Haisai^(-1)*eval(TiDu_cal)';
Haisai_cal='Haisai';
end
i=1;
while i < n
if abs(det(eval(Haisai_cal))) < 1e-6
disp('逆矩阵不存在！');
break;
end
x0=x0(：)+eval(dk);
if norm(eval(TiDu_cal)) < eps
X=x0;
result=eval(f_cal);
return;
end
i=i+1;
end
disp('无法收敛！');
X=[];
result=[];
end& r1 W: ? p+ |1 `* k4 {+ m

%% 外点惩罚函数-混合约束
syms x1 x2
f=(x1-2)^2+(x2-1)^2;
g=[-0.25*x1^2-x2^2+1];%修改成大于等于形式
h=[x1-2*x2+1];
x0=[2 2];
M=0.01;
C=3;
eps=1e-6;
[x,result]=waidian_hunhe(f,g,h,x0,M,C,eps,100)
function [x,result]=waidian_hunhe(f,g,h,x0,M,C,eps,k)
% f 目标函数
% g 不等式约束函数矩阵
% h 等式约束函数矩阵
% x0 初始值
% M 初始惩罚因子
% C 罚因子放大倍数
% eps 退出容差
% 循环次数
CF=sum(h.^2); %chengfa
n=1;
while n<k
%首先判断是不是在可行域内
gx=double(subs(g,symvar(g),x0));%计算当前点的约束函数值
index=find(gx<0);%寻找小于0的约束函数
F_NEQ=sum(g(index).^2);
F=matlabFunction(f+M*F_NEQ+M*CF);
x1=Min_Newton(F,x0,eps,100);
x1=x1'
if norm(x1-x0)<eps
x=x1;
result=double(subs(f,symvar(f),x));
break;
else
M=M*C;
x0=x1;
end
n=n+1;
end
end
3 `: C+ g* L1 \' `9 q# q1 O0 A! c

%% 内点惩罚函数
syms x1 x2 x
f=x1.^2+x2.^2;
g=[x1+x2-1;2*x1-x2-2];
x0=[3 1];
M=10;
C=0.5;
eps=1e-6;
[x,result]=neidian(f,g,x0,M,C,eps,100)
function [x,result]=neidian(f,g,x0,M,C,eps,k)
% f 目标函数
% g 不等式约束函数矩阵
% h 等式约束函数矩阵
% x0 初始值
% M 初始障碍因子
% C 障碍因子缩小倍数
% eps 退出容差
% k 循环次数
%惩罚项
Neq=sum((1./g));
n=1;
while n<k
F=matlabFunction(f+M*Neq);
[x1,result]=Min_Newton(F,x0,eps,100);
x1=reshape(x1,1,length(x0));
tol=double(subs(Neq,symvar(Neq),x1)*M);
if tol < eps
if norm(x1-x0) < eps
x=x1;
result=double(subs(f,symvar(f),x));
break;
else
x0=x1;
M=M*C;
end
else
if norm(x1-x0) < eps
x=x1;
result=double(subs(f,symvar(f),x));
break;
else
x0=x1;
M=M*C;
end
end
n=n+1;
end
end
. u# l9 k# E2 u2 Z

作者: yin123 时间: 2021-3-29 13:11
基于matlab约束优化之惩罚函数法

欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)