本站微信 收藏本站 设为首页
 找回密码
 注册
关于网站域名变更的通知
EDA365电子论坛网»电子社区 › 研发技术版块 › MATLAB论坛 › 基于matlab约束优化之惩罚函数法
返回列表 发新帖
查看: 570|回复: 1
打印 上一主题 下一主题

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

[复制链接]
mutougeda
  • TA的每日心情
    怒
    2019-11-20 15:22

    • 签到天数: 2 天

    [LV.1]初来乍到
    跳转到指定楼层
    1#
    发表于 2021-3-29 11:22 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式

    EDA365欢迎您登录！

    您需要 登录 才可以下载或查看，没有帐号？注册

    x

    ! i$ u" w2 I" a+ h) D( q一、算法原理7 o" G. C. ?( X5 M+ ]1 F
    1、问题引入/ C, U# F. b( D7 @4 J
    之前我们了解过的算法大部分都是无约束优化问题，其算法有：黄金分割法，牛顿法，拟牛顿法，共轭梯度法，单纯性法等。但在实际工程问题中，大多数优化问题都属于有约束优化问题。惩罚函数法就可以将约束优化问题转化为无约束优化问题，从而使用无约束优化算法。/ ]' j' z( k8 X. ]% e: c- r
    - j& P# q  P( p# ]
    2、约束优化问题的分类3 x% Y9 j) E1 P% _$ _" k* ~
    约束优化问题大致分为三类：等式约束、不等式约束、等式+不等式约束。
    $ K/ t2 m, I/ g8 s7 T) o/ {0 L2 E# A8 Z4 v  ~$ b$ M; v- [
    其数学模型为：8 X: k/ N2 q$ Q3 h# u6 ~
    , T6 {  ]' v( X0 K
    等式约束9 J" }# q$ c4 x- K: g2 ^0 s' s

    1 |4 ~4 Q0 j6 L! D: amin f(x),x\in R^n
    9 j% |1 }7 f* R/ h& G6 t
    ( f- d/ @4 h* `* }2 K% r; Ks.t hv(x)=0,v=1,2,…,p<n
    , r: O, b5 o" f1 s" S! u9 K* g+ J5 m2 T; ]$ g: F$ c3 \. I& ]
    不等式约束
    : U" G  E' a; b" T/ I+ R
    5 C: O) _% Y4 d9 o' B2 xmin f(x),x\in R^n
    7 g* p& E2 ?9 b8 M2 X" i; I
    * R8 n4 p: s" i9 S2 Xs.t gu(x)\leqslant 0,u=1,2,…,p<n
    % L, n  N! K# @0 Q3 K1 |/ K9 s- J0 V8 C" O% H/ {" z& j  Q
    等式+不等式约束问题
    , P0 J% o* d; y: g- F3 r0 t  \2 _4 W- P/ k2 f+ X  A. P. g0 `* ?
    min f(x),x\in R^n
    : L/ ~- W: i$ ~% f  n4 {% C# t3 M# u4 T4 g
    s.t hv(x)=0,v=1,2,…,p<n
    . A) f! l) ?% i* p
    * `8 Y  L- t0 Vs.t gu(x)\leqslant 0,u=1,2,…,p<n
    # e# m' T5 ^& g: }1 ]6 `8 f9 m9 i" B/ B0 ^: E9 B
    3、惩罚函数法定义
    + \( N& g  C+ [, P/ \6 m  y/ c惩罚函数法（SUMT法）又称序列无约束极小化技术，将等式约束与不等式约束的条件，经过适当定义的复合函数加到原目标函数上构造了惩罚函数，从而取消了约束，转而求解一系列无约束优化问题。; V# u! z& N' _5 R* l" f

    1 F9 e& i2 v4 \  C0 Y. e$ Z0 M按照惩罚函数再优化过程中的迭代点是否在约束条件的可行域内，又分为内点法、外点法和混合法
    $ g0 x7 @; A* Z. u# o# p& H: {% E8 f7 ?( K; ?# D0 k$ d1 }
    内点法：迭代点再约束条件的可行域之内，只用于不等式约束。
    # Q1 I: {; V9 p# ~& N' @7 E; o: v& ?: x
    ! n2 h# k1 N$ M外点法：迭代点再约束条件的可行域之外，既用于不等式约束又可用于等式约束。
    2 R' Y1 _. T. z4 v7 N
    ( e% t4 w  ]: e) ^" I) O% b7 G3 Q4、外点惩罚函数法* U" e9 Q! l" w$ }: @( ]3 Z% F. }
    等式约束：
    4 z. L; d6 x$ b, X
    % n8 P7 \) T# g" P( @min f=x1.2+x2.2,x\in R^n
    & y# C" T% m! @: v
    / a8 L. x- d' X7 P" \# H0 `s.t h1(x)=x1-2=0,h2(x)=x2+3=0
    ) ^# w) |9 x3 }' H. F* L
    , F9 S- p& ?" @; f+ O) w算法步骤
    ! `+ i/ H' X& C" [$ Z: {( U0 `; b9 K" J: [+ n; K& x( J
    a、构造惩罚函数：F=f+M * { [ h1(x) ]^2 + [ h2(x) ]^2 } ，式中M为初始惩罚因子；4 I1 f5 z4 N' Q: F$ W3 c0 @

    6 W/ e. J) h% W( g" S/ B* v' Mb、然后用无约束优化极值算法求解(牛顿法)；1 c8 |9 D; G& T3 R
    8 W+ A9 d5 }" I% Q1 z
    c、 如果相邻两次惩罚函数无约束最优点之间的距离足够小【norm(x1-x0)<eps】，则收敛；
    3 U* Q0 x" h$ X# o' x  V3 y  V2 d; M5 s" R
        否则放大惩罚因子M=C*M，式中C为 罚因子放大系数；
    + |$ _/ b3 z3 ~  {* d% c8 L( Y
    8 [7 y3 j6 z/ x6 I) D3 V& {  }
    • d、转步骤a继续迭代；
      % I1 ~3 {/ m3 G4 N8 J9 e

    , U6 W3 `0 y% C  d7 l4 [$ B, |0 B/ u+ Q

    3 s3 p9 c' `8 [4 o' R2 h! M二、源代码
    % A# `% @' \7 p' {
    • %% 外点惩罚函数法-等式约束
    • 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; {% N$ L8 M7 X6 J
       
    ' @$ b4 {- `  Z: Q9 }: O. l# `" c3 T9 r3 a

    $ B: ^2 e$ U, [& m' N1 D- G9 F3 F
    , r" _# s6 m1 e! l
    0 K1 ?" ~- L4 }* W) ?& K6 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
      + s. S" I# x- [- i+ G

    6 [+ K8 b: M5 d& I+ i) N
    . b7 j) r1 \. n    - P: r& e, q0 {( a* J" j% Z4 v

    1 U) P: h: F6 D) D9 o
    - d2 o$ c& x! l
    • %% 内点惩罚函数
    • 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
      , o5 C; K6 p) G, i" k8 d
       * `; }: e' o  n3 P5 E. T: ~0 M& N
    : Q9 L/ v3 Y* P
    / k% F9 r  x0 `8 m# A$ M) J

    下载资料威望不够？点击查看获取威望的N种方法>>
    回复

    举报

    EDA365微信号及QQ群号！
    yin123
  • TA的每日心情
    怒
    2019-11-29 15:37

    • 签到天数: 1 天

    [LV.1]初来乍到
    2#
    发表于 2021-3-29 13:11 | 只看该作者
    基于matlab约束优化之惩罚函数法

    下载资料威望不够？点击查看获取威望的N种方法>>
    回复 支持 反对

    举报

    返回列表 发新帖
    高级模式
    B Color Image Link Quote Code Smilies
    您需要登录后才可以回帖 登录 | 注册

    本版积分规则

    关闭

    推荐内容上一条 /1 下一条

    EDA365公众号

    关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

    GMT+8, 2025-11-23 21:31 , Processed in 0.156250 second(s), 26 queries , Gzip On.

    深圳市墨知创新科技有限公司

    地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

    快速回复 返回顶部 返回列表