一、算法原理

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

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

其数学模型为：

等式约束

s.t

不等式约束

s.t

等式+不等式约束问题

s.t

s.t

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

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

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

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

等式约束：

s.t

算法步骤

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

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

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

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

d、转步骤a继续迭代；

matlab代码

关于牛顿法的具体讲解，给出该博客网址https://blog.csdn.net/STM89C56/article/details/105643162

%% 外点惩罚函数法-等式约束

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)

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

不等式约束：

s.t

算法步骤

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

,(外电惩罚函数，迭代点再可行域之外，不等式约束才起作用)

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

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

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

d、转步骤a继续迭代；

matlab代码

牛顿法代码与等数约束相同这里不再给出

%% 外点惩罚函数法-不等式约束

syms x1 x2

f=x1.^2+(x2-2).^2;% x1-1>=0,x2-2>=0

g=[x1-1;-x2+2];%修改成大于等于形式

x0=[0 0];

M=0.03;

C=3;

eps=1e-6;

[x,result]=waidian_Neq(f,g,x0,M,C,eps,100)

function [x,result]=waidian_Neq(f,g,x0,M,C,eps,k)

% f 目标函数

% g 不等式约束函数矩阵

% x0 初始值

% M 初始惩罚因子

% C 罚因子放大倍数

% eps 退出容差

% k 循环次数

n=1;

while n

%首先判断是不是在可行域内

gx=double(subs(g,symvar(g),x0));%计算当前点的约束函数值

index=find(gx<0);%寻找小于0的约束函数

F_NEQ=sum(g(index).^2);

F=matlabFunction(f+M*F_NEQ);

x1=Min_Newton(F,x0,eps,100);

x1=reshape(x1,1,length(x0))

if norm(x1-x0)

x=x1;

result=double(subs(f,symvar(f),x));

break;

else

M=M*C;

x0=x1;

end

n=n+1;

end

混合约束：

s.t

算法步骤

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

,(外电惩罚函数，迭代点再可行域之外，不等式约束才起作用)

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

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

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

d、转步骤a继续迭代；

matlab代码

%% 外点惩罚函数-混合约束

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

%首先判断是不是在可行域内

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)

x=x1;

result=double(subs(f,symvar(f),x));

break;

else

M=M*C;

x0=x1;

end

n=n+1;

end

end

s.t

算法步骤

a、构造惩罚函数：F=f+M * 1/{ g1(x)  +  g2(x) } ，式中M为初始惩罚因子,

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

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

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

d、转步骤a继续迭代；

matlab代码

%% 内点惩罚函数

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

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