怎样去实现基于matlab约束优化的惩罚函数法呢

　　一、算法原理
　　1、问题引入
　　之前我们了解过的算法大部分都是无约束优化问题，其算法有：黄金分割法，牛顿法，拟牛顿法，共轭梯度法，单纯性法等。但在实际工程问题中，大多数优化问题都属于有约束优化问题。惩罚函数法就可以将约束优化问题转化为无约束优化问题，从而使用无约束优化算法。
　　2、约束优化问题的分类
　　约束优化问题大致分为三类：等式约束、不等式约束、等式+不等式约束。
　　其数学模型为：
　　等式约束
　　min f（x），xin R^n
　　s.t hv（x）=0，v=1，2，…，p《n
　　不等式约束
　　min f（x），xin R^n
　　s.t gu（x）leqslant 0，u=1，2，…，p《n
　　等式+不等式约束问题
　　min f（x），xin 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，xin 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继续迭代；
　　matlab代码
　　二、源代码
　　%% 外点惩罚函数法-等式约束
　　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
　　
　　%% 外点惩罚函数-混合约束
　　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
　　
　　%% 内点惩罚函数
　　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