【状态估计】非线性受控动力系统的线性预测器——Koopman模型预测MPC（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>非线性受控动力系统的线性预测器——Koopman模型预测MPC研究1. Koopman理论的基本概念与非线性系统线性化Koopman算子是一种通过无限维线性空间描述非线性动力系统的数学工具其核心思想是通过可观测函数的演化实现全局线性化。该理论包含以下关键要素Koopman本征函数作为基向量将状态投影到高维空间形成线性动力学的基础。Koopman特征值表征系统的稳定性与动态特性如振荡频率、衰减速率其幂次运算可直接预测多步状态演化。Koopman模式用于将高维投影空间的信息还原到原始状态空间完成非线性系统的重构。通过扩展动态模式分解EDMD或深度学习如深度神经网络可从数据中学习有限维Koopman算子的近似表示将非线性系统转化为高维线性模型。例如在四旋翼无人机控制中通过引入旋转矩阵作为可观测函数可在SE(3)流形上实现全局线性化避免了欧拉角奇异性问题。2. 模型预测控制MPC的基本原理与非线性挑战MPC通过滚动优化有限时域内的控制输入序列实现对系统约束和性能指标的综合优化。其核心步骤包括预测模型构建基于系统动态方程预测未来状态。在线优化求解以当前状态为初始条件求解最优控制序列通常为二次规划问题。反馈校正仅实施优化序列的第一个控制量并在下一周期重新求解。非线性MPC的瓶颈在于非凸优化问题导致计算复杂度高且难以保证全局最优解。模型失配风险复杂非线性系统的精确建模困难影响预测精度。3. Koopman-MPC框架的构建与优势将Koopman模型与MPC结合形成数据驱动的线性预测-优化架构其实现流程如下数据驱动的Koopman建模通过EDMD或神经网络从实验/仿真数据中提取Koopman算子的有限维近似。选择物理可观测函数如能量、旋转矩阵以增强模型泛化能力。线性MPC设计基于Koopman线性模型构建状态空间方程将非线性MPC转化为凸二次规划问题。约束条件如输入限幅、状态边界可直接映射到高维线性空间。实时控制与在线更新结合卡尔曼滤波或模糊补偿处理模型-真实动态的失配问题。在线自适应方法如目标网络软更新可动态调整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. 挑战与未来研究方向模型维度与泛化能力高维Koopman模型可能导致过拟合需结合稀疏性促进方法如L1正则化。混合建模物理模型数据驱动可提升小样本下的泛化性能。在线学习与鲁棒性非平稳动态下的模型漂移问题需引入增量学习或记忆回放机制。结合鲁棒MPC如Zonotopic管状控制处理不确定性和外部扰动。硬件部署优化边缘计算架构下的分布式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)、标准差 RMSE1sqrt(mean((Cy*x_true-C1lift(1,:)*x1lift).^2)); RMSE2sqrt(mean((Cy*x_true-C2lift(1,:)*x2lift).^2)); RMSE3sqrt(mean((Cy*x_true-C3lift(1,:)*x3lift).^2)); std_1std(x1lift); std_2std(x2lift); std_3std(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 AliftA1lift; BliftB1lift; liftFunliftFun1; xliftx1lift; if RMSE1 RMSE3 fprintf(选择薄板样条径向基函数) else AliftA3lift; BliftB3lift; liftFunliftFun3; xliftx3lift; fprintf(选择多项式基函数) end else AliftA2lift; BliftB2lift; liftFunliftFun2; xliftx2lift; if RMSE2 RMSE3 fprintf(选择高斯基函数) else AliftA3lift; BliftB3lift; liftFunliftFun3; xliftx3lift; fprintf(选择多项式基函数) end end %% ********************** Feedback control ******************************** %比较了基于Koopman算子的MPC控制器(K-MPC)和基于局部线性化的MPC控制器(L-MPC)在两种情况下的性能 disp(Press any key for feedback control) pause a0; %记录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*(nhm),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(i1); 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(i1)-0.4 yrr(i1)0.4) % if abs(yr - x_koop(2,1)) 0.01 %判别误差是否在允许范围 因为Cy [0 1]所以只需要比较x_koop的第二行的元素 % disp([当t,num2str(i.),s时误差超过0.005]) % a1; %记录一下表示有UKMPC % %使用从0到此时刻的操作数据,构建的新的Koopman模型 % X [ UX]; % Y [ UY]; % U 2*rand(1, i1) - 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,Nlift1: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 a1 %判断是否更新过KMPC % % Update Koopman MPC -- UKMPC % Upzetazeta; % 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(Nsim1,1),-k,linewidth,lw_koop-1); hold on p2 plot([0:Nsim]*deltaT,ymin*ones(Nsim1,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)开始optvalInf不可行 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到TmaxY从-0.6到0.7 set(gca,FontSize,20); xlabel(time [s]) ylabel(y) title(reference tracking.) % figure %图6 % p3 plot([0:Nsim]*deltaT,ones(Nsim1,1),-k,linewidth,lw_koop-1); hold on % p4 plot([0:Nsim]*deltaT,-ones(Nsim1,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(Nsim1,1),-k,linewidth,lw_koop-1); hold on p4 plot([0:Nsim]*deltaT,-ones(Nsim1,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代码、数据、文章