精华内容
下载资源
问答
  • MPC的基础原理并不复杂,在深入原理之前,借用Matlab MPC toolbox熟悉MPC控制架构。 2. MPC Designer 2.1 流程 以matlab cstr(continuous stirred-tank reactor)[1]为例。反应外界环境的温度Tc是输入控制量,输入...

    1.Introduction

    MPC的基础原理并不复杂,在深入原理之前,借用Matlab MPC toolbox熟悉MPC控制架构。

    1.1 MPC modeling

    MPC关于模型的定义如下:
    在这里插入图片描述
    用simulink 表示上面的关系:
    在这里插入图片描述
    根据matlab对MPC处理的介绍:

    The MPC controller performs all estimation and optimization calculations using a discrete-time, delay-free, state-space system with dimensionless input and output variables.

    MPCdelay\color{red}{MPC对输入输出的delay做了处理}:

    Delay removal — If the discrete-time model includes any input, output, or internal delays, the absorbDelay command replaces them with the appropriate number of poles at z = 0, increasing the total number of discrete states. The InputDelay, OutputDelay, and InternalDelay properties of the resulting state-space model are all zero.

    控制对象的模型(可以看成是开环模型):
    在这里插入图片描述
    输入扰动模型(闭环模型的一部分):
    在这里插入图片描述
    输出扰动模型(闭环模型的一部分)
    在这里插入图片描述
    测量噪声扰动模型(闭环模型的一部分)
    在这里插入图片描述

    1.2 MPC State Prediction

    MPC是闭环控制模型,控制器状态,不仅只考虑开环模型的状态变量,还需要考虑各种扰动和测量误差,用于预测系统的下一个状态。
    在这里插入图片描述
    MPC控制器用于后续计算的状态方程稍不同于控制对象的状态方程:
    在这里插入图片描述

    • 控制器的状态变量由下面的状态变量组成:
      在这里插入图片描述
    • 控制器的输入变量:
      在这里插入图片描述
    • 更新后的状态矩阵如下:
      在这里插入图片描述
    • 加上状态观测器
      在这里插入图片描述
      在这里插入图片描述

    1.3 quadratic program

    MPC将线性模型转换成二次优化问题,基本的组成元素有:
    在这里插入图片描述
    cost的组成
    在控制系统设计的时候,需要综合考虑各种指标,因而cost也是由多个项组合而成
    在这里插入图片描述

    • out reference tracking: 跟踪的精度
      在这里插入图片描述
      在这里插入图片描述
    • Manipulated Variable Tracking
      输入权重,系统希望控制器输出的控制量越小越好(耗能更小…),如LQR控制器中的cost中有uRuTuRu^T项。下面的公式中,往往uj,target=0u_{j,target}=0
      在这里插入图片描述
      在这里插入图片描述
    • Manipulated Variable Move Suppression
      很多系统中倾向于控制变量的改变要比较smooth,尽量不出现较大的波动
      在这里插入图片描述
      在这里插入图片描述
    • softening constraints
      softeningconstraints\color{red}{softening \, 能解决特殊情况下需要突破constraints的情况}
      在这里插入图片描述
      在这里插入图片描述

    完整的形式:
    在这里插入图片描述
    其中Vj,minVj,maxV_{j,min}和V_{j,max}是ECR value
    在这里插入图片描述
    QR的矩阵
    假设状态转移矩阵如下:
    在这里插入图片描述
    转移矩阵如下:
    在这里插入图片描述
    在prediction horizon上y的预测数值:
    在这里插入图片描述
    简化这个公式:
    在这里插入图片描述
    系数矩阵为:
    在这里插入图片描述
    公式中Δu\Delta_u需要转换成控制变量ziz_i
    在这里插入图片描述
    举例:如果控制变量ziΔuz_i和\Delta_u的关系如下图。
    在这里插入图片描述
    设置Δu(0)=z0,Δu(2)=z1,Δu(5)=z2\Delta u(0)=z0, \Delta u(2)=z1, \Delta u(5)=z2
    在这里插入图片描述
    将状态方程进行简化:
    在这里插入图片描述
    在这里插入图片描述
    化成二次优化的标准形式:
    在这里插入图片描述
    相关权重则表示为:
    在这里插入图片描述
    各种限制:
    在这里插入图片描述

    2. MPC Designer

    2.1 流程

    以matlab cstr(continuous stirred-tank reactor)[1]为例。反应外界环境的温度Tc是输入控制量,输入液体浓度CAiC_{Ai}可能有一定的波动(并且没有办法去测量);水箱中的温度T是可观测的输出量,CAC_A是无法观测的输出量;并且水箱温度T和CAC_A也是状态变量。
    在这里插入图片描述
    在这里插入图片描述
    用一阶线性微分方程近似这个模型:
    在这里插入图片描述

    • step1:构建状态空间模型并设置模型具体结构
      在这里插入图片描述
      设置具体结构:
      在这里插入图片描述
    • step2:定义每个状态的物理意义
      在这里插入图片描述
    • step3:设置限制和权重,hard constraints 和softening constraints 在二次优化方程中的形式如下。
      在这里插入图片描述
      在这里插入图片描述
    • step4: 设置目标函数中的权重系数
      注意输出变量CAC_A因为观测变量,所以权重为0.
      在这里插入图片描述
    • step5: 设置输入信号,对控制器进行调试
      • 1.输入阶跃信号
        在这里插入图片描述
      • 2.测试系统抗扰动的能力
        输入一个阶跃的输入扰动,系统很快适应下来。
        在这里插入图片描述
    • step6: 调试MPC控制器参数
      同样的根据[2],prediction horizon 定义process error的权重,control horizon 定义输入权重。
      本例中,control horizon 是 prediction horizon的1/51/5,在process weigh的设置的时候,给process error的权重0.3。
      在这里插入图片描述

    2.2 samples

    2.2.1 单输入单输出样例

    构建一个单输入单输出的二阶惯性模型:
    在这里插入图片描述
    在matlab里面,只要配置好相关的输入输出,模型和限制就可以。

    plant = tf(1,[1 0 0]);
    Ts = 0.1;   
    p = 10;
    m = 3;
    mpcobj = mpc(plant, Ts, p, m);
    mpcobj.MV = struct('Min',-1,'Max',5);   % 因为只有一个输入变量,只需要设置该变量的限制
    mdl = 'mpc_doubleint';
    open_system(mdl)
    sim(mdl)
    

    其他的内容:包括卡尔曼估计隐藏状态,以及matlab如何构建QP equation是看不到的。
    在这里插入图片描述

    对比PID控制器和MPC控制器的效果,单输入单输出模型,MPC的优势并不是太大。(下图中红色的是MPC,蓝色是PID)
    在这里插入图片描述
    比较MPC和PID控制器的输出指令,下图中(黄色的是PID 蓝色的是MPC)。
    在这里插入图片描述

    2.2.2 多输入单输出

    假设多输入多输出的状态转移矩阵为:
    在这里插入图片描述

    % 设置MPC模型
    sys = ss(tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}));
    Ts = 0.2;
    model = c2d(sys,Ts);   
    model = setmpcsignals(model,'MV',1,'MD',2,'UD',3);  % 3个输入变量,5个状态变量
    % 设置MPC控制器
    mpcobj = mpc(model,Ts,10,3);  % 控制器参数设置
    mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10);   % 设置输入限制
    % 设置其他通道
    mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);   % 噪声模型
    % 设置sim
    Tstop = 30;                                       % simulation time
    Tf = round(Tstop/Ts);                        % number of simulation steps
    r = ones(Tf,1);                                  % reference signal
    v = [zeros(Tf/3,1);ones(2*Tf/3,1)];   % measured disturbance signal,如何作用在系统上
    sim(mpcobj,Tf,r,v)
    % 让仿真的效果更接近实际
    SimOptions = mpcsimopt(mpcobj);   % 设置MPC的相关参数
    d = [zeros(2*Tf/3,1);-0.5*ones(Tf/3,1)];        
    SimOptions.Unmeas = d;                          % unmeasured input disturbance signal
    SimOptions.OutputNoise=.001*(rand(Tf,1)-.5);    % output measurement noise
    SimOptions.InputNoise=.05*(rand(Tf,1)-.5);      % noise on manipulated variables
    % 手动修改kalman filter 用于估计无法直接观测的变量
    [L,M,A1,Cm1] = getEstimator(mpcobj);
    e = eig(A1-A1*M*Cm1)
    poles = [.8 .75 .7 .85 .6 .81];
    L = place(A1',Cm1',poles)';
    M = A1\L;
    setEstimator(mpcobj,L,M);
    

    上面这个系统用simulink模拟如下图:
    在这里插入图片描述
    打印出模型的详细信息:
    在这里插入图片描述
    MDUD\color{red}{输入扰动MD和UD是如何作用在这个系统上的}

    2.2.3 非线性的多输入多输出问题

    传递模型是一个非线性模型,采用EKF的思路,将当前时刻下的状态线性处理:
    在这里插入图片描述
    在这里插入图片描述

    plant = linearize('mpc_nonlinmodel');
    plant.InputName = {'Mass Flow';'Heat Flow';'Pressure'};  
    plant.OutputName = {'Temperature';'Level'};
    plant.InputUnit = {'kg/s' 'J/s' 'Pa'};
    plant.OutputUnit = {'K' 'm'};
    Ts = 0.2;                          
    p = 5;
    m = 2;
    mpcobj = mpc(plant,Ts,p,m);
    mpcobj.MV = struct('Min',{-3;-2;-2},'Max',{3;2;2},'RateMin',{-1000;-1000;-1000});
    mpcobj.Weights = struct('MV',[0 0 0],'MVRate',[.1 .1 .1],'OV',[1 1]);
    outdistmodel = tf({1 0;0 1},{[1 0 0 0],1;1,[1 0 0 0]});
    setoutdist(mpcobj,'model',outdistmodel);
    mdl2 = 'mpc_nonlinear_setoutdist';
    open_system(mdl2)
    sim(mdl2)
    

    Reference

    [1] https://www.youtube.com/watch?v=J6RqQ9qDDeU
    [2] https://www.youtube.com/watch?v=gMOcBSmjdkQ

    展开全文
  • mpc 模型预测 matlab 代码
  • 1.Introduction MPC算法本身非常容易理解,关键在于如何将工程中的物理对象转换成合适的数学模型。 2.MPC Design 2.1时间参数 系统有控制 3.Applications 3.1倒立摆模型

    1.Introduction

    MPC算法本身非常容易理解,关键在于如何将工程中的物理对象转换成合适的数学模型。

    2.MPC Design

    2.1时间参数

    系统有控制周期TsT_s,prediction horizon pp和control horizon mm,这三个参数根据系统的开环属性进行设置,一般不会通过这三个参数调整系统特性。

    • 控制周期
      控制频率应至少要大于控制对象带宽的4倍,控制周期越短越增加了控制器的运算量。理论上周期周期越短,对扰动的抑制作用越强,如果将控制周期缩短一倍,但是对相同的扰动的抑制,没有明显的增加,就不需要减短。
    • prediction horizon
      如果期望的闭环响应时间是T,则T ≈ pTs
    • control horizon
      一般control horizon 比较小1-3

    2.2 constraints

    添加hard constraints 可能会造成infeasible solution

    • hard constraints
      MV 和MV Rate 的constraints 尽量不要全部设置成hard
      OV的constraints尽量全部设置成soft
    • soft constraints ECR

    0 — No violation allowed (hard constraint)
    0.05 — Very small violation allowed (nearly hard)
    0.2 — Small violation allowed (quite hard)
    1 — average softness
    5 — greater-than-average violation allowed (quite soft)
    20 — large violation allowed (very soft)

    • soft constraints weight
      比其他三项cost 之和要大1-2个数量级,如果violation 太大了,在这个基础上可以适当的扩大2-5倍。
    • time variant constraints 和 time variant weight
      对于时变的限制和时变的权重,一般给prediction中的每个点分别分配限制和权重,但是这样调试起来及其麻烦。

    2.3 调试权重

    2.3.1调试工具:

    • sensitivity:
    [J,sens] = sensitivity(MPCobj,PerfFunc,PerfWeights,Tstop,r,v,simopt,utarget)
    

    将主要项的cost记录下来,J是所有项的cost
    在这里插入图片描述
    在这里插入图片描述

    • review
      将MPC控制做一个综述报告
      在这里插入图片描述

    2.3.2 weight tuning

    matlabnnyc\color{red}{这个部分matlab还是有很多地方解释的非常不清楚:n是什么,nyc是什么}
    在这里插入图片描述

    2.4 plant with delay

    • 设置模型
      如果模型满足下面的关系
      在这里插入图片描述
    % 用传递模型表示:只通过简单的输入输出的关系也可以作为系统模型
    g11 = tf(12.8,[16.7 1],'IOdelay',1.0,'TimeUnit','minutes');
    g12 = tf(-18.9,[21.0 1],'IOdelay',3.0,'TimeUnit','minutes');
    g13 = tf(3.8,[14.9 1],'IOdelay',8.1,'TimeUnit','minutes');
    g21 = tf(6.6,[10.9 1],'IOdelay',7.0,'TimeUnit','minutes');
    g22 = tf(-19.4,[14.4 1],'IOdelay',3.0,'TimeUnit','minutes');
    g23 = tf(4.9,[13.2 1],'IOdelay',3.4,'TimeUnit','minutes');
    DC = [g11 g12 g13;
          g21 g22 g23];
    DC.InputName = {'Reflux Rate','Steam Rate','Feed Rate'};
    DC.OutputName = {'Distillate Purity','Bottoms Purity'};
    DC = setmpcsignals(DC,'MD',3);
    mpcDesigner(DC)
    

    关于调试参数
    PMtd,max/ΔtP-M \gg t_{d,max}/\Delta t
    td,maxt_{d,max}是最大延时时间;
    Δt\Delta t是控制周期

    2.5 测试MPC控制器的鲁棒性

    matlab的案例是比较,同一个控制器面对建模误差的鲁棒性。

    • 模型误差
    A = [-0.0285 -0.0014; -0.0371 -0.1476];
    B = [-0.0850 0.0238; 0.0802 0.4462];
    C = [0 1; 1 0];
    D = zeros(2,2);
    CSTR = ss(A,B,C,D);
    CSTR.InputName = {'T_c','C_A_i'};
    CSTR.OutputName = {'T','C_A'};
    CSTR.StateName = {'C_A','T'};
    CSTR = setmpcsignals(CSTR,'MV',1,'UD',2,'MO',1,'UO',2);
    
    dB = [0 0;0.05 0];
    perturbUp = CSTR;
    perturbUp.B = perturbUp.B + dB;
    
    perturbDown = CSTR;
    perturbDown.B = perturbDown.B - dB;
    
    • 观测step的反应
    step(CSTR,perturbUp,perturbDown)
    legend('CSTR','peturbUp','perturbDown')
    

    曲线的差别:
    在这里插入图片描述

    • 在mpcDesigner中添加噪声,比较同一个控制器对不同模型的抗噪声的能力
    • 在这里插入图片描述

    在这里插入图片描述

    3.Applications

    3.1 dc motor model

    在这里插入图片描述
    转换成物理模型:
    在这里插入图片描述
    matlab 上模型为:
    在这里插入图片描述
    模型状态变量是:
    在这里插入图片描述
    状态转移矩阵:
    在这里插入图片描述
    具体模型参数:输出参数[θ,T][\theta, T]θ\theta是measured outputs, TT是unmeasured outputs
    在这里插入图片描述
    最终的跟踪效果:
    在这里插入图片描述

    3.2 Paper Machine Process

    化工造纸的过程:
    在这里插入图片描述
    状态变量x:
    在这里插入图片描述
    在这里插入图片描述
    输入变量:[Gp,Gw,Np,Nw]{[G_p, G_w, N_p, N_w]}
    控制变量Gp,GwG_p, G_w
    在这里插入图片描述
    Np(consistency of stock entering the feed tank)是measured disturbance,
    Nw(the white water consistency)是unmeasured disturbance
    输出变量: [H2,N1,N2]{[H_2, N_1, N_2]}
    状态转移矩阵:
    在这里插入图片描述
    MPC模型:
    在这里插入图片描述
    跟踪的效果:
    在这里插入图片描述

    3.3倒立摆模型

    倒立摆是典型的非线性模型:

    • 线性模型的初始化:
    opspec = operspec(mdlPlant);
    opspec.States(1).Known = true;
    opspec.States(1).x = 0;
    opspec.States(3).Known = true;
    opspec.States(3).x = 0;
    
    • 将倒立摆模型进行线性化
    io(1) = linio([mdlPlant '/dF'],1,'openinput');
    io(2) = linio([mdlPlant '/F'],1,'openinput');
    io(3) = linio([mdlPlant '/Pendulum and Cart System'],1,'openoutput');
    io(4) = linio([mdlPlant '/Pendulum and Cart System'],3,'openoutput');
    options = findopOptions('DisplayReport',false);
    op = findop(mdlPlant,opspec,options);
    plant = linearize(mdlPlant,op,io);
    plant.InputName = {'dF';'F'};
    plant.OutputName = {'x';'theta'};
    pole(plant)
    plant = setmpcsignals(plant,'ud',1,'mv',2);
    

    mpc模型参数:
    在这里插入图片描述

    展开全文
  • matlab MPC 工具箱的英文介绍,很详细

    热门讨论 2012-05-26 14:53:49
    应该是官方的使用说明 很完整 具体
  • mpc MATLAB 代码

    2014-06-05 09:22:27
    mpc MATLAB 代码
  • MATLAB_MPC.zip

    2020-03-06 12:19:48
    基于状态空间模型的模型预测控制,搭建simulink平台,使用.m文件编辑S函数,实现多步模型预测控制。
  • matlab中的mpc工具箱

    2018-04-18 18:32:41
    matlab中的mpc工具箱
  • MPC工具箱提供的MPC模块不能实现权重参数的实时修改,有必要自己编写一个实现模型预测控制算法的matlab function。 主义事项 以下物理量必须设置相同,不然容易报错或求解不出理论控制量: 1. matlab func模块的...

    MPC工具箱提供的MPC模块不能实现权重参数的实时修改,有必要自己编写一个实现模型预测控制算法的matlab function。

    • 主义事项

    求解QP问题的时候使用哪一个函数更好????

    mpcqpsolver   (To be removed) Solve a quadratic programming problem using the KWIK algorithm

    quadprog

    mpcInteriorPointSolver

    以下物理量必须设置相同,不然容易报错或求解不出理论控制量:

    1. matlab func模块的调用频率

    2. MPC算法的采样频率。

    3. 被控对象的状态空间方程的离散化频率。

    以一个双积分系统为例子,MPC的实现代码如下

    function U = fcn(Ref,x,v,q1,ut)
    % clear
    % s=tf('s')
    % G=1/s^2
    % G=ss(G)
    % A=G.A
    % B=G.B
    % C=eye(2)
    % D=zeros(2,1)
    % sys=ss(A,B,C,D)
    % plant=c2d(sys,0.01)
    coder.extrinsic('quadprog');
    U = 0;
    deltaU =0;
    
    Ak = ...
      [1     0;
      0.01     1];
     
    Bk = ...
    [0.01;
    5e-05];
    
    Ck=eye(2);
    Dk=zeros(2,1);
    
    dk=zeros(6,1);
    ek=zeros(4,1);
    % 权重矩阵
    Q=diag([1.5+q1,100]);  % 1.52 100  2.8 
    R=diag([0]);
    S=diag([0]);
    pho=0.05;
    e=0.1;
    
    Ts=0.01;
    
    xk=[v;x;ut];
    yref=Ref;
    
    % 维数说明
    m=1; % 控制
    n=2; % 状态变量
    n0=m+n; % 新的状态维度
    Hc=2;
    Hp=10;
    p=size(Ck,1) ; %输出量维数
    % 已知量
    u_1 = ut;
    U_1=kron(ones(Hc,1),u_1);
    % syms deltau deltau1 deltau2 deltau3 deltau4
    % deltaU=[deltau;deltau1;deltau2;deltau3;deltau4]
    % U=M*deltaU+U_1
    
    %% 令 x=[x(k);u(k-1)]  控制量变为 delta u  则新系统的状态空间为
    Ak=[ Ak,Bk;zeros(m,n),eye(m) ];
    Bk=[Bk;eye(m)];
    Ck=[Ck,Dk];
    Dk=Dk;
    dk=[dk;zeros(m,1)];
    ek=ek;
    %% 预测输出矩阵 k+1 ... K+Hp    共Hp个
    % 控制输出矩阵  k...K+Hc-1   共Hc个
    
    % compute PSIk
    [m1,n1]=size(Ck*Ak);
    PSIk=zeros(m1*Hp , n1);
    for i=1:Hp
        PSIk( (i-1)*m1+1 :(i-1)*m1+m1 , 1:n1)=Ck*Ak^i;
    end
    % compute THETAk
    [m2,n2]=size(Ck*Bk);
    THETAk=zeros(m2*Hp , n2*Hc);
    for i=1:Hp
        for j=1:Hc
            if i>=j
                THETAk( (i-1)*m2+1:(i-1)*m2+m2, (j-1)*n2+1:(j-1)*n2+n2 ) = Ck*Ak^(i-j)*Bk;
            elseif j-i==1
                THETAk( (i-1)*m2+1:(i-1)*m2+m2, (j-1)*n2+1:(j-1)*n2+n2)=Dk;
            else
                THETAk( (i-1)*m2+1:(i-1)*m2+m2, (j-1)*n2+1:(j-1)*n2+n2)=zeros(m2,n2);
            end
        end
    end
    % compute TAUk
    [m3,n3]=size(Ck);
    TAUk=zeros(m3*Hp , n3*Hp);
    for i=1:Hp
        for j=1:Hp
            if i>=j
                TAUk( (i-1)*m3+1:(i-1)*m3+m3, (j-1)*n3+1:(j-1)*n3+n3 ) = Ck*Ak^(i-j);
            else
                TAUk( (i-1)*m3+1:(i-1)*m3+m3, (j-1)*n3+1:(j-1)*n3+n3)=zeros(m3,n3);
            end
        end
    end
    % compute PHIk
    PHIk=zeros(n0*Hp,1);  
    % compute LAMBDA
    LAMBDAk=zeros(p*Hp,1);   %% 线性时变时候不为0!!!!!!!!!!!!!
    
    % compute Qe
    [m4,n4]=size(Q);
    Qe=zeros(m4*Hp,n4*Hp);
    for i=1:Hp
        for j=1:Hp
            if i==j
                Qe( (i-1)*m4+1:(i-1)*m4+m4, (j-1)*n4+1:(j-1)*n4+n4 ) = Q;
            else
                Qe( (i-1)*m4+1:(i-1)*m4+m4, (j-1)*n4+1:(j-1)*n4+n4 )=zeros(m4,n4);
            end
        end
    end
    % compute Re
    [m5,n5]=size(R);
    Re=zeros(m5*Hc,n5*Hc);
    for i=1:Hc
        for j=1:Hc
            if i==j
                Re( (i-1)*m5+1:(i-1)*m5+m5, (j-1)*n5+1:(j-1)*n5+n5 ) = R;
            else
                Re( (i-1)*m5+1:(i-1)*m5+m5, (j-1)*n5+1:(j-1)*n5+n5 )=zeros(m5,n5);
            end
        end
    end
    % compute Se
    [m6,n6]=size(S);
    Se=zeros(m6*Hc,n6*Hc);
    for i=1:Hc
        for j=1:Hc
            if i==j
                Se( (i-1)*m6+1:(i-1)*m6+m6, (j-1)*n6+1:(j-1)*n6+n6 ) = S;
            else
                Se( (i-1)*m6+1:(i-1)*m6+m6, (j-1)*n6+1:(j-1)*n6+n6 )=zeros(m6,n6);
            end
        end
    end
    
    K=tril(ones(Hc));
    Im=eye(m,m);
    M=kron(K,Im);
    
    [m7,n7]=size(yref);
    Yrefk=zeros(m7*Hp,1);
    for i=1:Hp
        Yrefk( (i-1)*m7+1:(i-1)*m7+m7 , 1)=yref;
    end
    epsilon=PSIk*xk+TAUk*PHIk+LAMBDAk-Yrefk;
    %% 最终变量 Hk Gk Pk
    Hk=...
        [ 2*(THETAk'*Qe*THETAk+Re+M'*Se*M) ,  zeros(m*Hc,1) ;
        zeros(1,m*Hc)                    ,  pho             ];
    
    Gk=...
        [  2*epsilon'*Qe*THETAk+2*U_1'*Se*M  ,  0  ];
    
    
    Pk=...
        [  epsilon'*Qe*epsilon + U_1'*Se*U_1 + pho*e^2  ];
    
    
    %% QP问题
    H=Hk;
    c=Gk';
    A=[];
    b=[];
    Aeq=[];
    beq=[];
    Minv=M^-1;
    VLB=[Minv*(-ones(m*Hc,1)-U_1);-1];
    VUB=[Minv*(ones(m*Hc,1)-U_1);1];
    [x,z,fla,out]=quadprog(H,c,A,b,Aeq,beq,VLB,VUB);
    deltaU = x(1);
    
    U = ut(1)+deltaU  ;
    
    
    


    MPC算法的参考文献为:

    基于模型预测控制的无人驾驶车辆轨迹跟踪控制算法研究_孙银健

     

    !!!!!

    声明:若果您觉得您被侵犯了权利,请联系本博客作者,同样的,若引用了本文也请注明出处。

    !!!!!

    需要源代码、slx文件或技术支持的童鞋请联系我。

     

     

     

     

     

    展开全文
  • 自编Matlab代码实现MPC基本算法

    万次阅读 多人点赞 2019-12-28 06:27:53
    #用Matlab实现MPC的基本算法 写之前我在网上找了一下相关的文章,发现没有一个写得很直观的答案,刚好最近在研究MPC,于是决定自己写一个,力求做到直接了当,让入门的同学一眼就能看明白。 首先要介绍一下一个基于...

    自编Matlab代码实现MPC的基本算法

    本文属于入门级的教程,适合刚入坑的萌新看。在写之前我在网上找了一下相关的文章,发现没有一个写得很直观的答案,刚好最近在研究MPC,于是决定自己写一个,针对最简单的线性约束系统,推导其优化的过程原理,揭示其简单的结构。另外关于MPC的理论,可以自行参考相关教材,这里不再赘述。

    首先要介绍一下一个基于Matlab的工具MPT3,这是一个优化工具箱,同样也能调用MPC,感兴趣的同学可以自己点进去看一下。今天要用的模型也是来自这里面的一个demo,可以参考这个网址https://www.mpt3.org/UI/RegulationProblem

    一、模型准备

    这里用最简单的线性模型来进行解释,用差分方程进行描述如下:
    x+=Ax+Bu x^{+}=Ax+Bu
    其中,A=[1101]A=\begin{bmatrix}1& 1\\ 0& 1\end{bmatrix}B=[10.5]B=\begin{bmatrix} 1\\ 0.5\end{bmatrix}. 初始状态x0=[40]x_0=\begin{bmatrix} 4 &0\end{bmatrix}'. 然后是状态量以及控制量的约束限制:
    [55]x[55],1u1 \begin{bmatrix} -5\\ -5\end{bmatrix}\leq x\leq \begin{bmatrix} 5\\ 5\end{bmatrix}, -1\leq u\leq1
    一般的线性MPC中用的是二次优化,分别包含了针对状态量的代价函数xkQxkx_k'Qx_k项以及针对控制量的代价函数ukRuku_k'Ru_k:
    J(k,Np)=mint=kk+Np(xtQxt+utRut) J_{(k,N_p)}=\mathrm {min} \sum_{t=k}^{k+N_p}(x_t'Qx_t+u_t'Ru_t)
    表示的是在时刻tt,滚动时域步长为NpN_p的优化目标函数,目的是要让状态量xx以及控制量uu尽量的小,即使用最少的能量使系统状态量尽快稳定到零点(或者是稳定到某一个状态,这也是一般优化问题的目的)。其中代价系数矩阵的取值分别为:
    Q=[1001],R=1 Q=\begin{bmatrix}1& 0\\ 0& 1\end{bmatrix}, R=1
    其取值可以根据设计要求在状态量与控制量的大小之间平衡。预测时域的长度选为Np=5N_p=5,模型及参数准备完毕开始进行下一步。

    二、两个关键问题

    2.1 未来状态量的变量替换

    为了使用Matlab中的quadprog函数进行在线二次优化,求解Matlab中的优化问题,需要对优化目标进行转换,变换为关于控制量utu_t的优化函数,回顾优化目标函数,里面包含了在NpN_p步后的未来系统状态量xtx_t,这里需要将其描述为关于优化目标utu_t的函数。由于是线性系统,通过不停的迭代可以很容易得到xtx_t关于utu_t的表达式:
    xk+1k=Axk+Bukxk+2k=Axk+1k+Buk+1=A2xk+ABuk+Buk+1xk+Npk=Axk+Np1k+Buk+Np1=ANpxk+ANp1Buk++ABuk+Np2+Buk+Np1 x_{k+1|k}=Ax_{k}+Bu_k \\ x_{k+2|k}=Ax_{k+1|k}+Bu_{k+1}=A^2x_k+ABu_k+Bu_{k+1} \\ …… \\ x_{k+N_p|k}=Ax_{k+N_p-1|k}+Bu_{k+N_p-1}=A^{N_p}x_{k}+A^{N_p-1}Bu_k+……+ABu_{k+N_p-2}+Bu_{k+N_p-1}
    如果感兴趣可以自己尝试推导一下。然后将其描述为矩阵形式:
    [xk+1kxk+2k...xk+Npk]=[AA2...ANp]xk+[B0...0ABB...0............ANp1BANp2B...B][ukuk+1...uk+Np1] \begin{bmatrix}x_{k+1|k}\\ x_{k+2|k}\\ ...\\ x_{k+N_p|k}\end{bmatrix}=\begin{bmatrix}A\\ A^2\\ ...\\ A^{N_p}\end{bmatrix}x_k+\begin{bmatrix}B& 0 & ... &0 \\ AB& B & ... &0 \\ ...& ... & ... & ...\\ A^{N_p-1}B& A^{N_p-2}B & ... &B \end{bmatrix}\begin{bmatrix}u_k\\ u_{k+1}\\ ...\\ u_{k+N_p-1}\end{bmatrix}
    简化表示为:
    X=A~xk+B~U X=\tilde{A}x_k+\tilde{B}U
    其中X=[xk+1k,xk+2k,...,xk+Npk],U=[uk,uk+1,...,uk+Np1]X=[x'_{k+1\mid k},x'_{k+2\mid k},...,x'_{k+N_{\textrm{p}}\mid k}]', U=[u'_k,u'_{k+1},...,u'_{k+N_{\textrm{p}}-1}]'A~,B~\tilde{A},\tilde{B}为对应的矩阵,这样就将将来的状态量表示为了关于优化目标UU的表达式,可以代入优化目标中进行求解了。这里值得注意的是在Matlab中,通过适当的技巧可以得到At,BtA_t,B_t,具体方法可以参考后文的代码。

    2.2 优化目标函数的转换

    在完成了未来状态量的变量替换之后,进一步将优化目标函数替换为全部关于优化目标utu_t以及当前状态量xkx_k的表达式。将前文中的优化目标函数重新表示为:
    J(k,Np)=min XQ~X+UR~U J_{(k,N_p)}=\mathrm {min} \ X'\tilde{Q}X+U'\tilde{R}U
    其中Q~,R~\tilde{Q},\tilde{R}是适当增广变换之后的优化代价系数矩阵,感兴趣的同学可以自己推导。具体的变化方法在后文Matlab代码中可见。总之就是尽量将优化目标简化,使其能够直接传递到quadprog函数的参数中。参考Matlab中quadprog函数的形式,调用方法为:
    [x,fval,exitflag]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)minx 12xHx+fxs.t.    Axb                  Aeqx=beq                 lbxub [x,fval,\mathrm {exitflag}]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x_0,options) \\ \underset{x}{\mathrm {min}} \ \frac{1}{2}x'Hx+f'x \\ s.t. \ \ \ \ Ax\leq b \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Aeq*x=beq \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ lb\leq x \leq ub
    为了将优化目标函数表示为标准形式,将xtx_t的表达式代入其中得到:
    J(k,Np)=(A~xk+B~U)Qt(A~xk+B~U)+URtU=U(B~QtB~+Rt)U+2xkA~Q~B~U+xkA~Q~A~xk J_{(k,N_p)}=(\tilde{A}x_k+\tilde{B}U)'Q_t(\tilde{A}x_k+\tilde{B}U)+U'R_tU \\ =U'(\tilde{B}'Q_t\tilde{B}+R_t)U+2x'_k\tilde{A}'\tilde{Q}\tilde{B}U+x_k'\tilde{A}'\tilde{Q}\tilde{A}x_k
    其中,UU为优化目标,第三项与优化目标无关,故可以将其省略掉,所以有
    H=2(B~Q~B~+R~)f=(2xkA~Q~B~) H=2(\tilde{B}'\tilde{Q}\tilde{B}+\tilde{R}) \\ f=(2x_k'\tilde{A}'\tilde{Q}\tilde{B})
    再根据状态量xtx_t的约束得到A,b(注意这里的A要与系统系数矩阵区别开来)的表达式:
    5A~xk+B~U5{B~U5A~xkB~U5+A~xk[B~B~]U[5A~xk5+A~xk] -5\leq \tilde{A}x_k+\tilde{B}U\leq 5 \\ \Rightarrow \left\{\begin{matrix} \tilde{B}U\leq 5-\tilde{A}x_k\\ -\tilde{B}U\leq 5+\tilde{A}x_k \end{matrix}\right. \\ \Rightarrow \begin{bmatrix} \tilde{B}\\ -\tilde{B}\end{bmatrix}U\leq \begin{bmatrix}5-\tilde{A}x_k\\5+\tilde{A}x_k\end{bmatrix}
    故可得到:
    A=[B~B~],b=[5A~xk5+A~xk] A=\begin{bmatrix} \tilde{B}\\ -\tilde{B}\end{bmatrix},b=\begin{bmatrix}5-\tilde{A}x_k\\5+\tilde{A}x_k\end{bmatrix}
    同样地,可以由输入的约束轻易得到lb,ublb,ub。至此准备工作完成,可以开始用Matlab实现MPC的优化过程了。

    三、Matlab在线优化及代码

    clear; clc; close all;
    % 线性系统系数矩阵
    A=[1 1; 0 1]; B=[1; 0.5];
    % 初始状态量-如果不能在下一步回到约束范围内,则会造成无解
    x0=[4; 0];
    % 预测步长
    Np=5;
    % 优化目标参数,加权矩阵
    Q=eye(2); R=1;
    % 转化为用控制量ut表示的,关于状态量的推导方程的矩阵
    At=[]; Bt=[]; temp=[];
    % 转换后的加权矩阵
    Qt=[]; Rt=[];
    
    % 加权矩阵的计算过程,以及推导方程矩阵的叠加过程
    for i=1:Np
            At=[At; A^i];
            Bt=[Bt zeros(size(Bt,1), size(B,2));
                A^(i-1)*B temp];
            temp=[A^(i-1)*B temp];
            
            Qt=[Qt zeros(size(Qt,1),size(Q,1));
                zeros(size(Q,1),size(Qt,1)) Q];
            Rt=[Rt zeros(size(Rt,1),size(R,1));
                zeros(size(R,1),size(Rt,1)) R];
    end
    % 控制量ut的上下限
    lb=-1*ones(Np,1);
    ub=1*ones(Np,1);
    % 控制量ut的初始值
    u0=zeros(Np,1);
    % 转换后的优化目标函数矩阵,循环优化函数中H后的表达式为优化目标的另一项
    H=2*(Bt'*Qt*Bt + Rt);
    % 转换后的优化中的不等式约束左边系数矩阵,后面循环中的bi为不等式右边
    Ai=[Bt; -Bt];
    % 声明u来保存每一步采用的控制量
    u=[];
    x=x0;
    xk=x0;
    
    
    for k=1:50
        
        % 关于ut的不等式约束,实际上约束的是状态量,常数4就是状态量约束的上下边界
        bi=[5-At*xk; 5+At*xk];
        % 一切准备就绪,进行二次优化
        [ut, fval, exitflag]=quadprog(H,(2*xk'*At'*Qt*Bt)',Ai,bi,[],[],lb,ub,u0);
        % 显示求解结果是否正常
        fprintf('%d\n', exitflag)
    
        % 采用优化得到的控制量的第一个元素作为实际作用的控制量,代入到原系统中得到下一个时刻的状态量
        u(k) = ut(1);
        x(:, k+1) = A*x(:, k) + B*u(k);
        xk = x(:, k+1);
        % 对优化初始值进行修改,采用预测值的后段作为下一步的初始值
        u0 = [ut(2:Np); ut(Np)];
        
    end
    
    
    figure();
    plot(x');
    legend('x_1','x_2');
    
    figure();
    plot(u);
    legend('u');
    
    

    以上为全部的代码,以下为仿真得到的结果:
    在这里插入图片描述
    在这里插入图片描述

    四、总结

    本文只是针对线性系统的MPC做了仿真,并且简化了很多内容。在后期大家可以根据自身的需求,在此基础上进行修改,比如考虑非线性系统,增加输出约束,增加终端代价函数和终端约束来保证稳定性等。或者进一步考虑跟踪,改变预测步长NpN_p来观察其对优化性能的影响,改变优化目标中的代价矩阵观察控制量的变化等等。尽情探索吧,如果有时间我会再更新进一步关于MPC的内容。最后欢迎大家来一起讨论!

    更新

    应一些朋友的要求,在本文的基础上又更新了一篇关于定点跟踪的文章——自编Matlab代码实现MPC定点跟踪,基本原理不变,对本文的代码做了一些修改,增加了一个作图的功能,感兴趣的朋友可以看看。

    展开全文
  • matlab开发-MPC自动变速器MPCInsimulinkv2。在Simulink中使用MPC的教程。
  • MPC模型预测控制(二)-MATLAB代码实现

    万次阅读 多人点赞 2018-12-24 15:18:34
    update:MPC的QQ群 https://blog.csdn.net/tingfenghanlei/article/details/85046120在这篇文章里主要讲了下MPC的原理和C++实现的一个简单例子。 这篇文章里主要写MPCMATLAB实现...我看MPCMATLAB代码实现...
  • MATLAB环境中,实现基于动态矩阵控制(DMC)的模型预测控制(MPC)的源程序代码。 在MATLAB环境中,实现基于动态矩阵控制(DMC)的模型预测控制(MPC)的源程序代码。
  • matlab开发-MPC555电机控制功能块组。MPC555目标的附加I/O块-专门针对TPU功能
  • 自编Matlab代码实现MPC定点跟踪

    千次阅读 2020-05-15 20:10:28
    自编Matlab代码实现MPC(线性)定点跟踪 这篇文章是在上一篇的基础上进一步对自编的MPC-Matlab代码进行了一定程度的延伸,主要是增加了定点跟踪功能,以及将每一步的预测状态量与控制量作图的功能。 一、模型准备及...
  • matlab code for mpc controller. 有step response 实例和 经典水箱加热 范例
  • 请问大家,在使用matlabmpc工具箱的非线性模型预测控制中nlmpc函数时 ,如果系统矩阵A中含有随时间变化的已知参数该如何处理呢?希望有知道的朋友解答解答。跪谢。
  • 离散控制Matlab代码自主车辆MPC C ++中的自行车模型上的MPC 控制方法: 模型预测控制=后视最优控制+约束最优控制 在以下三种情况下,使用模型预测控制器以恒定的前进速度实现了车辆的横向控制: 直线控制 换道操纵 ...
  • matlab运行神经网络的代码低数据限制下用于模型预测控制的非线性动力学的稀疏识别 带有控制的非线性动力学的稀疏识别(SINDYc)与模型预测控制(MPC)相结合。 该框架通过少量测量学习受外源控制变量影响的非线性...
  • MPC工具箱MATLAB .pdf

    2020-03-27 11:15:04
    simulink中MPC工具箱的简单介绍,英文版的,pdf文档,英文版虽然看起来有点费劲,但还是有用的,里面有实例。

空空如也

空空如也

1 2 3 4 5 ... 18
收藏数 348
精华内容 139
关键字:

matlabmpc

matlab 订阅