💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
目录
💥1 概述
📚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
📚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 参考文献
部分理论来源于网络，如有侵权请联系删除。
