【状态估计】非线性受控动力系统的线性预测器——Koopman模型预测MPC(Matlab代码实现)

 👨‍🎓个人主页

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势:🌞🌞🌞博客内容尽量做到思维缜密,逻辑清晰,为了方便读者。

⛳️座右铭:行百里者,半于九十。

📋📋📋本文目录如下:🎁🎁🎁

目录

💥1 概述

非线性受控动力系统的线性预测器——Koopman模型预测MPC研究

1. Koopman理论的基本概念与非线性系统线性化

2. 模型预测控制(MPC)的基本原理与非线性挑战

3. Koopman-MPC框架的构建与优势

4. 典型应用案例与性能对比

5. 挑战与未来研究方向

6. 结论

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、数据、文章


💥1 概述

文献来源:

本文提出了一类非线性受控动力系统的线性预测器。其基本思想是将非线性动力学提升(或嵌入)到其演化近似线性的高维空间中。在不受控制的情况下,这个过程相当于与非线性动力学相关的库普曼算子的数值近似。在这项工作中,我们将Koopman算子扩展到受控动力系统中,并应用扩展动态模态分解(EDMD)来计算算子的有限维近似,从而使该近似具有线性受控动力系统的形式。在数值实例中,以这种方式获得的线性预测器表现出优于现有线性预测器的性能,例如基于局部线性化或所谓的Carleman线性化的线性预测器。重要的是,构建这些线性预测器的过程完全是数据驱动的,而且非常简单——它归结为数据的非线性变换(提升)和提升空间中的线性最小二乘问题,可以很容易地解决大型数据集。这些线性预测器可以很容易地使用线性控制器设计方法来设计非线性动态系统的控制器。我们特别关注模型预测控制(MPC),并表明以这种方式设计的MPC控制器具有与具有相同数量的控制输入和相同维数的状态空间的线性动力系统的MPC相当的底层优化问题的计算复杂度。重要的是,在提出的MPC方案中,状态和控制输入的线性不等式约束以及状态的非线性约束可以以线性方式施加。类似地,状态变量中的非线性成本函数可以用线性方式处理。我们处理全状态测量情况和输入输出情况,以及具有干扰/噪声的系统。数值算例验证了该方法

原文摘要:

This paper presents a class of linear predictors for nonlinear controlled dynamical systems. The basicidea is to lift (or embed) the nonlinear dynamics into a higher dimensional space where its evolution is approximately linear. In an uncontrolled setting, this procedure amounts to numerical approximations of the Koopman operator associated to the nonlinear dynamics. In this work, we extend the Koopman operator to controlled dynamical systems and apply the Extended Dynamic Mode Decomposition (EDMD) to compute a finite-dimensional approximation of the operator in such a way that this approximation has the form of a linear controlled dynamical system. In numerical examples, the linear predictors obtained in this way exhibit a performance superior to existing linear predictors such as those based on local linearization or the so called Carleman linearization. Importantly, the procedure to construct these linear predictors is completely data-driven and extremely simple – it boils down to a nonlinear transformation of the data (the lifting) and a linear least squares problem in the lifted space that can be readily solved for large data sets. These linear predictors can be readily used to design controllers for the nonlinear dynamical system using linear controller design methodologies. We focus in particular on model predictive control (MPC) and show that MPC controllers designed in this way enjoy computational complexity of the underlying optimization problem comparable to that of MPC for a linear dynamical
system with the same number of control inputs and the same dimension of the state-space. Importantly, linear inequality constraints on the state and control inputs as well as nonlinear constraints on the state can be imposed in a linear fashion in the proposed MPC scheme. Similarly, cost functions nonlinear in the state variable can be handled in a linear fashion. We treat both the full-state measurement case and the input–output case, as well as systems with disturbances/noise. Numerical examples demonstrate the approach.1

非线性受控动力系统的线性预测器——Koopman模型预测MPC研究

1. Koopman理论的基本概念与非线性系统线性化

Koopman算子是一种通过无限维线性空间描述非线性动力系统的数学工具,其核心思想是通过可观测函数的演化实现全局线性化。该理论包含以下关键要素:

  • Koopman本征函数:作为基向量,将状态投影到高维空间,形成线性动力学的基础。
  • Koopman特征值:表征系统的稳定性与动态特性(如振荡频率、衰减速率),其幂次运算可直接预测多步状态演化。
  • Koopman模式:用于将高维投影空间的信息还原到原始状态空间,完成非线性系统的重构。

通过 扩展动态模式分解(EDMD) 或深度学习(如深度神经网络),可从数据中学习有限维Koopman算子的近似表示,将非线性系统转化为高维线性模型。例如,在四旋翼无人机控制中,通过引入旋转矩阵作为可观测函数,可在SE(3)流形上实现全局线性化,避免了欧拉角奇异性问题。

2. 模型预测控制(MPC)的基本原理与非线性挑战

MPC通过滚动优化有限时域内的控制输入序列,实现对系统约束和性能指标的综合优化。其核心步骤包括:

  • 预测模型构建:基于系统动态方程预测未来状态。
  • 在线优化求解:以当前状态为初始条件,求解最优控制序列(通常为二次规划问题)。
  • 反馈校正:仅实施优化序列的第一个控制量,并在下一周期重新求解。

非线性MPC的瓶颈在于:

  • 非凸优化问题:导致计算复杂度高且难以保证全局最优解。
  • 模型失配风险:复杂非线性系统的精确建模困难,影响预测精度。
3. Koopman-MPC框架的构建与优势

将Koopman模型与MPC结合,形成数据驱动的线性预测-优化架构,其实现流程如下:

  1. 数据驱动的Koopman建模
    • 通过EDMD或神经网络从实验/仿真数据中提取Koopman算子的有限维近似。
    • 选择物理可观测函数(如能量、旋转矩阵)以增强模型泛化能力。
  2. 线性MPC设计
    • 基于Koopman线性模型构建状态空间方程,将非线性MPC转化为凸二次规划问题。
    • 约束条件(如输入限幅、状态边界)可直接映射到高维线性空间。
  3. 实时控制与在线更新
    • 结合卡尔曼滤波或模糊补偿处理模型-真实动态的失配问题。
    • 在线自适应方法(如目标网络软更新)可动态调整Koopman模型参数。

核心优势

  • 计算效率提升:线性化模型使优化问题复杂度从O(n³)降为O(n²),实验表明控制器可运行于100Hz以上。
  • 全局稳定性保障:通过Lyapunov函数或鲁棒管状MPC设计,确保闭环系统的稳定性。
  • 多场景适应性:在软体机器人、自动驾驶、化工过程等复杂系统中验证了有效性。
4. 典型应用案例与性能对比
  • 软体机器人控制:基于Koopman-MPC的气动机械臂轨迹跟踪误差(1.26 cm)显著低于传统线性模型(2.45 cm)。
  • 四旋翼无人机:在SE(3)流形上实现8字形高速轨迹跟踪,计算速度达100Hz,误差低于非线性MPC。
  • 锂离子电池状态估计:深度Koopman卡尔曼滤波较传统方法估计精度提升30%以上。
  • 燃料电池压力控制:结合模糊补偿的Koopman-MPC在质子交换膜电解槽中有效抑制压力波动。
5. 挑战与未来研究方向
  1. 模型维度与泛化能力
    • 高维Koopman模型可能导致过拟合,需结合稀疏性促进方法(如L1正则化)。
    • 混合建模(物理模型+数据驱动)可提升小样本下的泛化性能。
  2. 在线学习与鲁棒性
    • 非平稳动态下的模型漂移问题需引入增量学习或记忆回放机制。
    • 结合鲁棒MPC(如Zonotopic管状控制)处理不确定性和外部扰动。
  3. 硬件部署优化
    • 边缘计算架构下的分布式Koopman-MPC实现。
    • FPGA/GPU加速的实时求解器开发。
6. 结论

Koopman-MPC通过非线性系统的全局线性化,成功将复杂控制问题转化为高效凸优化问题,在计算效率、控制精度和鲁棒性方面展现出显著优势。随着深度学习与自适应算法的进一步融合,其在智能制造、自动驾驶、能源系统等领域的应用前景广阔,但需在模型可解释性、实时学习机制和硬件集成层面持续突破。

📚2 运行结果

 

 

部分代码:

figure %图1
stairs([0:Nsim-1]*deltaT,u_dt(0:Nsim-1),'linewidth',2); hold on
title('Control input'); xlabel('time [s]')

figure %图2
lw_koop = 3;
plot([0:Nsim]*deltaT,Cy*x_true,'-b','linewidth', lw_koop); hold on
plot([0:Nsim]*deltaT,C1lift(1,:)*x1lift, '--g','linewidth',lw_koop) %thinplate
plot([0:Nsim]*deltaT,C2lift(1,:)*x2lift, ' -r','linewidth',lw_koop) %gauss
plot([0:Nsim]*deltaT,C3lift(1,:)*x3lift, '--c','linewidth',lw_koop) %polyharmonic
plot([0:Nsim]*deltaT,Cy*xloc, '--k','linewidth',lw_koop-1)
LEG = legend('True','thinplate','gauss','polyharmonic','Local at $x_0$');
set(LEG,'Interpreter','latex','location','n  ortheast','fontsize',30)
set(gca,'FontSize',25);
axis([0 1 -1.3 0.5])
xlabel('time [s]')
ylabel('y')
title('prediction comparison')

figure %图3
lw_koop = 3;
plot([0:Nsim]*deltaT,Cy*x_true - C1lift(1,:)*x1lift, '--g','linewidth',lw_koop) %thinplate
hold on
plot([0:Nsim]*deltaT,Cy*x_true - C2lift(1,:)*x2lift, ' -r','linewidth',lw_koop) %gauss
plot([0:Nsim]*deltaT,Cy*x_true - C3lift(1,:)*x3lift, '--c','linewidth',lw_koop) %polyharmonic
plot([0:Nsim]*deltaT,Cy*x_true - Cy*xloc, '--k','linewidth',lw_koop-1)
LEG = legend('thinplate','gauss','polyharmonic','Local at $x_0$');
set(LEG,'Interpreter','latex','location','n  ortheast','fontsize',30)
set(gca,'FontSize',25);
axis([0 1 -0.5 1])
xlabel('time [s]')
ylabel('y error')
title('prediction error comparison')


%4个基函数的Koopman模型预测比较的均方根误差(RMSE)、标准差
RMSE1=sqrt(mean((Cy*x_true-C1lift(1,:)*x1lift).^2));                      
RMSE2=sqrt(mean((Cy*x_true-C2lift(1,:)*x2lift).^2));  
RMSE3=sqrt(mean((Cy*x_true-C3lift(1,:)*x3lift).^2));  

std_1=std(x1lift);
std_2=std(x2lift);
std_3=std(x3lift);
sd=[std_1 std_2 std_3];

figure%图4
RMSE = categorical({'thinplate','gauss','polyharmonic'});
prices = [RMSE1 RMSE2 RMSE3];
bar(RMSE,prices);
xlabel('basic function')
ylabel('RMSE')
title('Average RMSE for prediction comparison')
%errorbar(prices,sd);

if RMSE1 <= RMSE2
    Alift=A1lift;
    Blift=B1lift;
    liftFun=liftFun1;
    xlift=x1lift;
    if  RMSE1 <= RMSE3
        fprintf('选择薄板样条径向基函数')
    else
      Alift=A3lift;
      Blift=B3lift;
      liftFun=liftFun3;
      xlift=x3lift;
      fprintf('选择多项式基函数')
    end
else 
    Alift=A2lift;
    Blift=B2lift;
    liftFun=liftFun2;
    xlift=x2lift;
   if RMSE2 <= RMSE3
      fprintf('选择高斯基函数')
   else
      Alift=A3lift;
      Blift=B3lift;  
      liftFun=liftFun3;
      xlift=x3lift;
      fprintf('选择多项式基函数')
   end
end


%% ********************** Feedback control ********************************
%比较了基于Koopman算子的MPC控制器(K-MPC)和基于局部线性化的MPC控制器(L-MPC)在两种情况下的性能
disp('Press any key for feedback control')
pause

a=0; %记录UKMPC
Tmax = 3; % Simlation legth
Nsim = Tmax/deltaT;
REF = 'cos'; % 'step' or 'cos'
switch REF % 参考输出,即可以测量到的数据
    case 'step' 
        ymin = -0.6; %constraint
        ymax = 0.6;
        x0 = [0;0.6];
        yrr = 0.3*( -1 + 2*([1:Nsim] > Nsim/2)  ); % 分段常数reference,no state constraints
    case 'cos'
        ymin = -0.4; %constraint
        ymax = 0.4;
        x0 = [-0.1;0.1];
        yrr = 0.5*cos(2*pi*[1:Nsim] / Nsim); % time-varying reference,constraints on the output imposed
end

% Define Koopman controller
C = zeros(1,Nlift); C(1) = 1;  %C如何选取
% Weight matrices 权重矩阵
Q = 1;
R = 0.01;
% Prediction horizon 预测范围
Tpred = 1;
Np = round(Tpred / deltaT); 
% Constraints 
xlift_min = [ymin ; nan(Nlift-1,1)];
xlift_max = [ymax ; nan(Nlift-1,1)];


% Build Koopman MPC controller -- KMPC
koopmanMPC = getMPC(Alift,Blift,C,0,Q,R,Q,Np,-1, 1, xlift_min, xlift_max,'qpoases');%----------这个函数里面包含了求解成本函数(目标函数)

% Initial condition for the delay-embedded state (assuming zero control in the past) 延迟嵌入状态的初始条件(假设过去的控制为零)
x = x0;
zeta0 = [Cy*x ; NaN(nd*(nh+m),1)];
for i = 1:nd
    uprev = zeros(m,1);
    xp = f_ud(0,x,uprev);
    zeta0 = [Cy*xp ; uprev ; zeta0(1:end-nh-m)];  %式子(9.22)
    x = xp;
end
x0 = x;
x_koop = x0; x_loc = x0;
zeta = zeta0; % Delay-embedded "state" 需要提升的状态

XX_koop = x0; UU_koop = [];
XX_loc = x0; UU_loc = [];
ZETA = [];
UX = [];UY = [];

% Get Jacobian of the true dynamics (for local linearization MPC) -- LMPC,得到真实动力学的雅可比矩阵(对于局部线性化MPC)
x = sym('x',[2 1]); syms u;
f_ud_sym = f_ud(0,x,u);
u_loc = 0;
Jx = jacobian(f_ud_sym,x);
Ju = jacobian(f_ud_sym,u);

wasinfeas= 0;
ind_inf = [];


% Closed-loop Simulation Start ------------ to get an optimal solution u.----就是Algorithm 1 :Lifting MPC – closed-loop operation
UU_koop1 = [];
for i = 0:Nsim-1
    if(mod(i,10) == 0)
        fprintf('Closed-loop simulation: iterate %i out of %i \n', i, Nsim)
    end
    
    % Current value of the reference signal 参考信号的值
    yr = yrr(i+1);
    YR = [yr];
    
    % Koopman MPC -- KMPC 式子(9.22)
    %UX = [UX zeta];
        z0 = liftFun(zeta); % Lift  算法第2步
        u_koop = koopmanMPC(z0,yr); % Get control input,the optimal solution u*   算法第3步
        x_koop = f_ud(0,x_koop,u_koop); % Update true state, apply u* to the system  算法第5步
        zeta = [ Cy*x_koop ; u_koop; zeta(1:end-nh-m)]; % Update delay-embedded state  在每个时间步骤中,从提升状态ψ(xk)初始化预测。
    %UY = [UY zeta];
%    
%     %unpdate Koopman predictor -- UKMPC
%     %如果系统的实际输出与参考输出之间的误差超过允许范围,则需要使用新的运行数据重新构建模型,更新koopmanMPC,更新A、B
%     if (yrr(i+1)>=-0.4 & yrr(i+1)<=0.4)
%         if abs(yr - x_koop(2,1)) >= 0.01  %判别误差是否在允许范围     因为Cy = [0 1],所以只需要比较x_koop的第二行的元素
%             disp(['当t=',num2str(i.'),'s时,误差超过0.005'])
%             a=1; %记录一下,表示有UKMPC
%            %使用从0到此时刻的操作数据,构建的新的Koopman模型
%             X = [ UX];
%             Y = [ UY];
%             U = 2*rand(1, i+1) - 1;       
%             UpXX_koop = [XX_koop]; %UKMPC之前的数据就继承KMPC的数据就行,这样画图也好画
%             UpUU_koop = [UU_koop];
%             %lift
%             UXlift = liftFun(X);
%             UYlift = liftFun(Y);
%             %find A、B
%             W = [UYlift ; X];
%             w = [UXlift ; U];
%             G = w*w';
%             V = W*w';
%             H = V * pinv(G);%pinv是求伪逆矩阵  % Matrix [A B;C 0]
%             Aulift = H(1:Nlift,1:Nlift);
%             Bulift = H(1:Nlift,Nlift+1:end);
%             %Update Koopman MPC controller
%             U_koopmanMPC = getMPC(Alift,Blift,C,0,Q,R,Q,Np,-1, 1, xlift_min, xlift_max,'qpoases');
%             fprintf('系统的实际输出与参考输出之间的误差超过允许范围,使用新的运行数据重新构建模型,更新koopmanMPC')
%         else
%         end
%     else
%     end
%     
%     %UKMPC
%     if a==1  %判断是否更新过KMPC
%         % Update Koopman MPC -- UKMPC
%         Upzeta=zeta;
%         Upx_koop = x_koop;
%         xulift = liftFun(zeta); % Lift
%         Upu_koop = U_koopmanMPC(xulift,yr); % Get control input
%         Upx_koop = f_ud(0,Upx_koop,Upu_koop); % Update true state
%         Upzeta = [Cy*Upx_koop ; Upu_koop; zeta(1:end-nh-m)];
%         % Store values
%         UpXX_koop = [UpXX_koop Upx_koop];
%         UpUU_koop = [UpUU_koop Upu_koop];
%     else
%     end
%     
    % Local linearization MPC -- LMPC
    Aloc = double(subs(Jx,[x;u],[x_loc;u_loc])); % Get local linearization
    Bloc = double(subs(Ju,[x;u],[x_loc;u_loc]));
    cloc = double(subs(f_ud_sym,[x;u],[x_loc;u_loc])) - Aloc*x_loc - Bloc*u_loc;
    [U_loc,~,optval] = solveMPCprob(Aloc,Bloc,Cy,cloc,Q,R,Q,Np,-1, 1,[nan;ymin],[nan;ymax],x_loc,yr); % Get control input
    u_loc = U_loc(1:m,1);
    if(optval == Inf) % Detect infeasibility
        ind_inf = [ind_inf i];disp(['当i=',num2str(i.'),'s时,该optval不可行'])
        wasinfeas = 1;
    end
    x_loc = f_ud(0,x_loc,u_loc); % Update true state
    
    % adaptive MPC -- AMPC 通过将每个控制区间的动力学依次线性化,建立预测模型
    % 还没想好怎么加这个对比
    
    % Store values
    XX_koop = [XX_koop x_koop]; 
    UU_koop = [UU_koop u_koop];
    XX_loc = [XX_loc x_loc];
    UU_loc = [UU_loc u_loc];
    ZETA = [zeta0 zeta];
    
end

if(isempty(ind_inf))
    ind_inf = Nsim;
end


%% Plot
% Control signal
figure %图5
% Output (y = x_2)
p1 = plot([0:Nsim]*deltaT,ymax*ones(Nsim+1,1),'-k','linewidth',lw_koop-1); hold on
p2 = plot([0:Nsim]*deltaT,ymin*ones(Nsim+1,1),'-k','linewidth',lw_koop-1);
p3 = plot([0:Nsim]*deltaT,XX_koop(2,:),'-b','linewidth',lw_koop); hold on %'KMPC'
p4 = plot([0:ind_inf(1)-1]*deltaT,XX_loc(2,1:ind_inf(1)),'--g','linewidth',lw_koop); %画的图,除去了不可行的部分(从ind_inf(1)开始,optval=Inf,不可行)
p5 = plot([0:Nsim-1]*deltaT,yrr,'--r','linewidth',lw_koop); %'Reference'
% p6 = plot([0:Nsim]*deltaT,UpXX_koop(2,:),':c','linewidth',lw_koop); hold on %'UKMPC'
LEG  = legend([p3,p4,p5,p2],'KMPC','LMPC','Reference','Constraints');
set(LEG,'Interpreter','latex','location','southeast')
set(LEG,'Fontsize',18)
axis([0,Tmax,-0.6,0.7]) %x从0到Tmax,Y从-0.6到0.7
set(gca,'FontSize',20);
xlabel('time [s]')
ylabel('y')
title('reference tracking.')

% figure %图6
% p3 = plot([0:Nsim]*deltaT,ones(Nsim+1,1),'-k','linewidth',lw_koop-1); hold on
% p4 = plot([0:Nsim]*deltaT,-ones(Nsim+1,1),'-k','linewidth',lw_koop-1); hold on
% p1 = plot([0:ind_inf(1)-1]*deltaT,UU_loc(1:ind_inf(1)),'--g','linewidth',lw_koop); hold on
% p2 = plot([0:Nsim-1]*deltaT,UU_koop,'-b','linewidth',lw_koop); hold on
% axis([0,Tmax,min(UU_loc) - abs(min(UU_loc))*0.1, max(UU_loc)*1.1] )
% LEG  = legend([p2,p1,p3],'K-MPC','L-MPC','Constraint');
% set(LEG,'Interpreter','latex','location','southeast')
% set(gca,'FontSize',31);
% xlabel('time [s]')
% ylabel('tracking error')

figure %图6 Control inputs
p3 = plot([0:Nsim]*deltaT,ones(Nsim+1,1),'-k','linewidth',lw_koop-1); hold on
p4 = plot([0:Nsim]*deltaT,-ones(Nsim+1,1),'-k','linewidth',lw_koop-1); hold on
p1 = plot([0:ind_inf(1)-1]*deltaT,UU_loc(1:ind_inf(1)),'--g','linewidth',lw_koop); hold on
p2 = plot([0:Nsim-1]*deltaT,UU_koop,'-b','linewidth',lw_koop); hold on
%p5 = plot([0:Nsim-1]*deltaT,UpUU_koop,':c','linewidth',lw_koop); hold on %'UKMPC' 图线和KMPC一样了,说明有问题
axis([0,Tmax,min(UU_loc) - abs(min(UU_loc))*0.1, max(UU_loc)*1.1] )
LEG  = legend([p2,p1,p3],'KMPC','LMPC','Constraint');
set(LEG,'Interpreter','latex','location','southeast')
set(gca,'FontSize',31);
xlabel('time [s]')
ylabel('u')
title('Control inputs')

🎉3 参考文献

部分理论来源于网络,如有侵权请联系删除。

🌈4 Matlab代码、数据、文章

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值