精华内容
下载资源
问答
  • Matlab约束优化

    2021-02-15 09:34:15
    约束优化

    约束优化方法

    约束优化是寻求具有约束条件的线性或者非线性规划问题解的数值算法。

    随机方向法

    基本思路是在可行域内选择一个起始点,产生若干个随机方向,从中选择一个可以使目标函数值下降最快的随机方向作为可行搜索方向d\vec{d}.
    从初始点x0x_0出发,沿着d\vec{d}方向以一定步长搜索,得到新点xx,并且满足约束条件
    {gj(x)0(j=1,2,,m)f(x)<f(x0) \begin{cases} g_j(x)\leq 0(j=1,2,\dots, m)\\ f(x)<f(x_0) \end{cases}

    复合形法

    nn维空间的可行域中选取KK个点作为初始复合形的顶点。然后比较复合形中各顶点目标函数的大小,其中目标函数最大的点为坏点,用其他点(映射点)替换坏点,形成新的复合形。

    可行方向法

    从可行点出发,沿着下降可行方向进行搜索,求出使目标函数值下降的新的可行点。

    二次规划(QP)

    QP数学模型为
    minx12xTHx+fTxs.t.{AxbAeqx=beqlbxub \min_x\frac{1}{2}x^THx+f^Tx\\ s.t. \begin{cases} Ax\leq b\\ Aeq\cdot x=beq\\ lb\leq x\leq ub \end{cases}
    求解指令
    [x, fval, exitflag, output, lambda]=quadprog(H, f, A, b, lb, ub, x0, options):返回解x和目标值fval,和解x处包含的拉格朗日乘子的lambda结构参数。

    拉格朗日法

    拉格朗日乘子法就是求函数f(x1,x2,)f(x_1, x_2, \dots)g(x1,x2,)=0g(x_1, x_2, \dots)=0的约束条件下的极值的方法。

    • PH算法(等式约束)
      效用函数为
      minM(x,u(k),σk)=f(x)+u(k)Th(x)+σkh(x)Th(x) \min M(x, u^{(k)}, \sigma_k)=f(x)+u^{(k)T}h(x)+\sigma_k h(x)^Th(x)
      判定函数为
      ϕk=h(x(k)) \phi_k=\lVert h(x^{(k)})\rVert
    • PHR算法
      松弛变量法
      M(u,v,p)=f(x)+12ρ2ρl{[max(0,ui+ρgi(x))]2ui2}+j=1lvjhj(x)+ρ2j=1lhj2(x) M(u, v, p)=f(x)+\frac{1}{2\rho}\sum_{2\rho}^l\{[\max(0, u_i+\rho g_i(x))]^2-u_i^2\}+\sum_{j=1}^l v_jh_j(x)+\frac{\rho}{2}\sum_{j=1}^l h_j^2(x)
      乘子的修正公式为
      {vj(k+1)=vj(k)+ρhj(x(k)),j=1,,lui(k+1)=max[0,ui(k)+ρgi(x(k))],i=1,,m \begin{cases} v_j^{(k+1)}=v_j^{(k)}+\rho h_j(x^{(k)}), j=1, \dots, l\\ u_i^{(k+1)}=\max[0, u_i^{(k)}+\rho g_i(x^{(k)})], i=1,\dots, m \end{cases}
      判断函数为
      ϕk={j=1lhj2(x(k))+i=1mmax(gi(x(k)),ui(k)ρ)2}1/2 \phi_k=\{\sum_{j=1}^l h_j^2(x^{(k)})+\sum_{i=1}^m\max(-g_i(x^{(k)}), \frac{u_i^{(k)}}{\rho})^2\}^{1/2}

    应用:投资组合管理

    %%
    clear all;
    clc;
    ExpReturn = [0.0017, 0.0015, 0.0005];
    ExpCovariance = [
        0.0010, 0.0004, 0.0005;
        0.0004, 0.0016, 0.0003;
        0.0005, 0.0003, 0.0013;
    ];
    
    NumPorts = 5;
    NumAssets = 3;
    PVal = 1;
    AssetMin = 0;
    AssetMax = [0.6, 0.7, 0.5];
    GroupA = [1 0 0];
    GroupB = [0 1 1];
    GroupMax = [0.7, 0.5];
    AtoBmax = 3;
    ConSet = portcons('PortValue', PVal, NumAssets, 'AssetLims', AssetMin, AssetMax, NumAssets, 'GroupComparison', GroupA, NaN, AtoBmax,...
        GroupB, GroupMax);
    
    %%
    A = ConSet(:, 1:end-1);
    b = ConSet(:, end);
    p = Portfolio;
    p = setAssetMoments(p, ExpReturn, ExpCovariance);
    p = setInequality(p, A, b);
    
    PortWts = estimateFrontier(p, NumPorts);
    [PortRisk, PortReturn] = estimatePortMoments(p, PortWts);
    
    disp('Risk');
    PortRisk
    disp('Return');
    PortReturn
    disp('Efficient Matrix');
    PortWts
    
    展开全文
  • MATLAB约束优化之惩罚函数法

    千次阅读 2020-04-25 11:14:21
    惩罚函数法就可以将约束优化问题转化为无约束优化问题,从而使用无约束优化算法。 2、约束优化问题的分类 约束优化问题大致分为三类:等式约束、不等式约束、等式+不等式约束。 其数学模型为: 等式约束...

    一、算法原理

    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继续迭代;

    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)<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

    不等式约束:

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

    s.t    g1(x)=x1-1\geqslant =0,h2(x)=x2-2\geqslant =0

    算法步骤

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

    u=1 ,g(x)<0;u=0,g(x)>=0,(外电惩罚函数,迭代点再可行域之外,不等式约束才起作用)

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

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

            否则放大惩罚因子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<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);
        x1=Min_Newton(F,x0,eps,100);
        x1=reshape(x1,1,length(x0))
        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

    混合约束:

    min f=(x1-2).^2+(x2-1).^2,x\in R^n

    s.t    g1(x)=-0.25*x1.^2-x2.^2+1\geqslant 0

    h2(x)=x1-2*x2+1 =0

    算法步骤

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

    u=1 ,g(x)<0;u=0,g(x)>=0,(外电惩罚函数,迭代点再可行域之外,不等式约束才起作用)

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

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

            否则放大惩罚因子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<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

    5、内点惩罚函数法

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

    s.t    g1(x)=x1+x2-1\geqslant 0

    g2(x)=2*x1-x2-2 \geqslant 0

    算法步骤

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

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

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

            否则缩小惩罚因子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<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

     

    展开全文
  • 惩罚函数法就可以将约束优化问题转化为无约束优化问题,从而使用无约束优化算法。 2、约束优化问题的分类 约束优化问题大致分为三类:等式约束、不等式约束、等式+不等式约束。 其数学模型为: 等式约束 s.t ...

    一、算法原理

    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继续迭代;

    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)<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

    不等式约束:

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

    s.t    g1(x)=x1-1\geqslant =0,h2(x)=x2-2\geqslant =0

    算法步骤

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

    u=1 ,g(x)<0;u=0,g(x)>=0,(外电惩罚函数,迭代点再可行域之外,不等式约束才起作用)

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

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

            否则放大惩罚因子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<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);
        x1=Min_Newton(F,x0,eps,100);
        x1=reshape(x1,1,length(x0))
        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

    混合约束:

    min f=(x1-2).^2+(x2-1).^2,x\in R^n

    s.t    g1(x)=-0.25*x1.^2-x2.^2+1\geqslant 0

    h2(x)=x1-2*x2+1 =0

    算法步骤

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

    u=1 ,g(x)<0;u=0,g(x)>=0,(外电惩罚函数,迭代点再可行域之外,不等式约束才起作用)

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

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

            否则放大惩罚因子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<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

    5、内点惩罚函数法

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

    s.t    g1(x)=x1+x2-1\geqslant 0

    g2(x)=2*x1-x2-2 \geqslant 0

    算法步骤

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

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

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

            否则缩小惩罚因子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<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

     

    展开全文
  • 一、算法原理 1、问题引入 之前我们了解过的算法大部分都是无约束优化问题,其算法有:黄金分割法,牛顿法,拟牛顿法,共轭梯度法,单纯性法等。但在实际工程问题中,大多数优化问题都属于有约束优化问题。惩罚函数...

    一、算法原理

    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继续迭代;

    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)<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
    

    四、备注

    完整代码或者代写添加QQ2449341593
    往期回顾>>>>>>
    【优化求解】基于matlab粒子群优化灰狼算法【含Matlab源码 006期】
    【优化求解】基于matlab多目标灰狼优化算法MOGWO 【含Matlab源码 007期】
    【优化求解】基于matlab粒子群算法的充电站最优布局【含Matlab源码 012期】
    【优化求解】基于matlab遗传算法的多旅行商问题【含Matlab源码 016期】
    【优化求解】基于matlab遗传算法求最短路径【含Matlab源码 023期】
    【优化求解】基于matlab遗传和模拟退火的三维装箱问题【含Matlab源码 031期】
    【优化求解】基于matlab遗传算法求解车辆发车间隔优化问题【含Matlab源码 132期】
    【优化求解】磷虾群算法【含matlab源码 133期】
    【优化求解】差分进化算法【含Matlab源码 134期】

    展开全文
  • Matlab约束优化

    2021-02-14 10:08:39
    约束优化及算法
  • MATLAB约束优化(UOM)

    2020-12-09 21:30:54
    MATLAB约束优化(UOM) 文章目录MATLAB约束优化(UOM)一、基本思想二、基本算法1、最速下降法(共轭梯度法)2、牛顿法3、拟牛顿法3.1 BFGS3.2 DFP三、MATLAB优化1、求解优化问题的主要函数2、优化函数的输入变量3、...
  • matlab开发-约束优化模拟退火。使用模拟退火进行连续约束优化
  • matlab开发-约束优化的CrowSearch算法。约束CSA
  • 约束优化matlab代码

    2018-04-19 17:16:46
    压缩包里有着关于无约束优化的代码,是利用matlab进行实现。
  • matlab 源代码 约束优化问题 经典奉献
  • 例4-1 MATLAB实现用M函数文件形式求解 syms s t; f=s^2+3*t^2-2*t*s+4*t+3*s; [x,minf]=minZBLH(f,[-2 6],[0.2 0.2],1.5,[t s],0.0001,0.0001) 坐标轮换minZBLH函数文件如下 function [x,minf] = minconPS2(f,x0,...
  • matlab实现约束优化——并行计算

    千次阅读 2017-12-06 20:16:39
    之前写完的约束优化代码需要在十八个测试问题上跑完,由于随机性的影响,需要进行多次测试问题,耗时太久,所以需要用到matlab的并行计算的功能。
  • 必须参数funx1x2 说明 fun 为目标函数用M文件或Inline定义 x1x2为目标函数的自变量的取值范围 options是一个结构型的变量用于指定优化参数可通过optimset函数设置 help optimset;例子1求解;例子2求解模型;方法1;2...
  • matlab解决无约束优化问题

    千次阅读 2020-06-06 21:01:21
    约束优化问题 要用到的数学知识: 1、向量范数与矩阵范数 2、多元函数梯度与Hessian阵 3、凸集与凸函数 特别要提示的是:如果该函数为凸函数,那么它有且仅有一个最优点,如果它的值不在无穷处,我们利用大部分...
  • MATLAB进行无约束优化

    2018-08-27 18:45:00
    首先先给出三个例子引入fminbnd和fminuc函数求解无约束优化,对这些函数有个初步的了解 求f=2exp(-x)sin(x)在(0,8)上的最大、最小值。 例2 边长3m的正方形铁板,四角减去相等正方形,制成方形无盖水槽。...
  • MATLAB求解优化问题的主要函数 类 模型 基本函数名 元函数极小|MinF(x)s.t.x1,x,x 无约束极小 Min F(X) X-fminunc(F, Xo) X-fminsearch'F 线性规划 Winc X X-linprog(c, A, b) s t AX-b 次规划 当nxx+c-x X=...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 805
精华内容 322
关键字:

matlab约束优化

matlab 订阅