基于matlab约束优化之惩罚函数法

mutougeda

" y6 \" a3 z2 o$ S' n& l3 i一、算法原理
6 u3 X( J! R3 P3 N  \: h" E. ?: ^3 c, d1、问题引入
之前我们了解过的算法大部分都是无约束优化问题，其算法有：黄金分割法，牛顿法，拟牛顿法，共轭梯度法，单纯性法等。但在实际工程问题中，大多数优化问题都属于有约束优化问题。惩罚函数法就可以将约束优化问题转化为无约束优化问题，从而使用无约束优化算法。

4 o" H: o! ]2 W& d- P1 H2、约束优化问题的分类
约束优化问题大致分为三类：等式约束、不等式约束、等式+不等式约束。

其数学模型为：

等式约束
5 |+ o$ W% d( Y3 ^2 ^6 D$ E( d
$ U2 z5 t  y/ |9 j4 J; x1 F) q& Imin f(x),x\in R^n
: C; O7 K. ?1 e0 E
s.t hv(x)=0,v=1,2,…,p<n

不等式约束

2 o6 H  T/ R  Z+ ^/ R: L% smin f(x),x\in R^n

( r6 I$ g: z2 y4 T% K! ms.t gu(x)\leqslant 0,u=1,2,…,p<n

. ~- E% O8 b/ p& s$ z等式+不等式约束问题
5 t6 I  j$ t1 q& Q5 W
& O$ R$ G; B5 H' h( Z2 Amin f(x),x\in R^n

3 I& Q/ ^5 n8 ?* L) x/ Ks.t hv(x)=0,v=1,2,…,p<n

( D, Q& Q# q* k5 ^' c( s: Us.t gu(x)\leqslant 0,u=1,2,…,p<n
- A/ S9 O6 {$ D) z
* N  k! }% ~+ z* |3、惩罚函数法定义
- [( k3 G0 T& v2 W惩罚函数法（SUMT法）又称序列无约束极小化技术，将等式约束与不等式约束的条件，经过适当定义的复合函数加到原目标函数上构造了惩罚函数，从而取消了约束，转而求解一系列无约束优化问题。
2 b& K0 `" Y) g6 s" \
+ s8 K. H* x' \$ Q% a; X按照惩罚函数再优化过程中的迭代点是否在约束条件的可行域内，又分为内点法、外点法和混合法
) O  F- ~6 x, d8 u
0 Y, n+ s! _8 D- K内点法：迭代点再约束条件的可行域之内，只用于不等式约束。

  t  j8 L9 D% O8 \( C7 Q外点法：迭代点再约束条件的可行域之外，既用于不等式约束又可用于等式约束。

4、外点惩罚函数法
等式约束：
4 K! l- ]5 L" C3 t
0 k4 h7 Q% P) s# i7 Omin f=x1.2+x2.2,x\in R^n

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

! S7 u1 H( W" n7 \( R, J# p# e( G/ f" n算法步骤

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

5 }: w) w# I* Db、然后用无约束优化极值算法求解(牛顿法)；
0 n& e# X# M, ~# T1 Y3 ^' V# B
( |. x& y5 D, S, A" H* [( L9 d- Oc、 如果相邻两次惩罚函数无约束最优点之间的距离足够小【norm(x1-x0)<eps】，则收敛；
" U- u* w* R) I5 l5 K5 p* I
3 a/ h: h, f; w( t, u    否则放大惩罚因子M=C*M，式中C为 罚因子放大系数；

• d、转步骤a继续迭代；
  8 ^0 p; P, U+ q4 b3 |5 n6 B

. D# I) s, d  A4 X' k$ B

0 u! P& @# T& s$ e5 [二、源代码
4 I% v. n) w, m  r3 e* Z

• %% 外点惩罚函数法-等式约束
• 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
  

   
9 N0 o7 J- d. E  i9 P
7 [0 C6 h) _* r6 {0 A) l


2 h! B& a1 G9 T9 `- f7 V& `

• %% 外点惩罚函数-混合约束
• 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
  

   
6 f. K2 T$ Z1 L: b

• %% 内点惩罚函数
• 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
  8 \5 S. @- q1 g! V

   
$ F+ y& c( p: y* d6 v. t
" x0 c: t( {) W, y  {2 ?, b" `

yin123

基于matlab约束优化之惩罚函数法