波浪能转换装置的模型预测控制策略研究（matlab代码实现）-CSDN博客

本文链接：https://blog.csdn.net/weixin_46039719/article/details/129803105

👨‍🎓个人主页：研学社的博客

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

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

⛳️座右铭：行百里者，半于九十。

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

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

💥1 概述

本文模拟从波浪能转换器（WEC）中提取的能量，当受控移动窗口阻塞 MPC 时，单设备。它还比较了使用标准MPC和GPC控制时WEC提取的能量。

摘要：海浪能是可再生能源最集中的来源之一。然而，到目前为止，它还没有达到商业化所需的经济可行性。为了提高波浪能转换器的效率，已经提出了几种先进的控制策略，包括模型预测控制（MPC）。然而，每个优化问题的计算负担都是传统（全自由度）MPC的缺点，这通常会限制其在系统实时控制中的应用。本文提出了一种移动窗口阻塞（MWB）方法，通过减少决策变量的数量来加快每个优化问题所需的时间。该方案控制的单器件点吸收器波能转换器的数值仿真证实了该方法的潜力。

📚2 运行结果

部分代码：

%% Clear All
clear all
close all
clc

%% Plot results from wecSim
% NEED file: wecSim_Output.mat
load('../data/WEC_ModelData.mat')
figure;
subplot(4,1,1)
plot(wave.time,wave.elevation)
hold on
plot(buoy.time,buoy.position)
title('Wave Elevation wecSim')
xlabel('time [s]')
ylabel('Elevation [m]')
grid on
grid minor
box on
subplot(4,1,2)
plot(buoy.time,buoy.position)
title('Displacement wecSim')
xlabel('time [s]')
ylabel('displacement [m]')
grid on
grid minor
box on
subplot(4,1,3)
plot(buoy.time,buoy.velocity)
title('Velocity wecSim')
xlabel('time [s]')
ylabel('velocity [m/s]')
grid on
grid minor
box on
subplot(4,1,4)
plot(buoy.time,buoy.ExcForce)
title('Force Excitation wecSim')
xlabel('time [s]')
ylabel('excForce [N]')
grid on
grid minor
box on

%% Build the Continous Matrices Ac, Bc and Cc for the Anaylisis of the WEC System
rho=1000; %[kg/m3] Water density
g= 9.81; %[m/s2] Acceleration of gravity
kh=rho*g*Sphere.wpArea; %[kg/s2] Hydrostatic Stiffness Coefficient
M=1 /(Sphere.Mass + Mu); %Sphere Mass
Ts=mean(diff(wave.time)); %Mean Sampling Time of Wave Force

%System matrix Ac in continuous time
Ac= [0 1 zeros(1,length(B_ss));
-M*kh 0 -M*C_ss;
zeros(length(B_ss),1) B_ss A_ss]; %Continous State Transition Matrix

%Input matrices Bc and Fc in continuous time
Bc = [0 M zeros(1,length(B_ss))]'; %Continous Input To State Matrix for the PTO force (Same for Wave Excitation Force)

%Output matrix Cc in continuous time
Cc=[1 0 zeros(1,length(B_ss));
0 1 zeros(1,length(B_ss))]; %Continous Output Matrix

%Form State Space
sysC = ss(Ac,Bc,Cc,0); %Declares Continuous State Space

%Discretize using ZOH
sysD = c2d(sysC,Ts,'zoh'); %Discretizes Using Zero-Order Hold (ZOH)

%% WEC Model
%System
Ad = sysD.A; %Discrete State Transition Matrix
Bd = sysD.B; %Discrete Input To State Matrix
Cd = sysD.C; %Discrete Output Matrix
Dd = sysD.D; %Discreete Input To Output Matrix
nx=size(Ad,1); %State Size
nu=size(Bd,2); %Input Size
ny=size(Cd,1); %Output Size

%% Delta MPC Setting
%MPC Parameters
Np=100; %Prediction Horizon
scaler=1000000; %Scaler for Matrix B

%Augmented system
A=[Ad,Bd*scaler,Bd*0;
Bd'*0,eye(nu),zeros(nu);
Bd'*0,zeros(nu),zeros(nu)]; %Augmented State Transition Matrix
Bw=[Bd*scaler;
zeros(nu);
zeros(nu)]; %Augmented Wave Input Matrix
B=[Bd*scaler;
eye(nu);
eye(nu)]; %Augmented Input Matrix
nxa=size(A,1); %Augmented State Size
C=[1,zeros(1,nxa-1);
0,1,zeros(1,nxa-2);
zeros(1,nxa-2),1,0;
zeros(1,nxa-1),1]; %Augmented Output Matrix

%Crossed Weights
Q=[0,0,0,0;
0,0,1,0;
0,1,0,0;
0,0,0,0]; %Output Peenalizations

%Constraints
DeltaUmax=2000000/scaler; %Input Increment Constraints (Scaled to avoid numerical problems)
DeltaUmin=-DeltaUmax;
Umax=200000/scaler; %Input Constraints (Scaled to avoid numerical problems)
Umin=-Umax;

%Matrices
[H,f1,f2,G,F,Fw]=MPC(A,B,Bw,C,Q,0,Np); %MPC Matrix Generation Function
H=(H'+H)/2; %Symmetrize Main Hessian

%% Lambda Penalized MPC Setting
lambda=1.12; %Input Penalization Term
Q_lambda=[0,0,0,0;
0,0,1,0;
0,1,lambda,0;
0,0,0,0]; %Output Penalizations w/Input Penalization
[H_lam,f1_lam,f2_lam]=MPC(A,B,Bw,C,Q_lambda,0,Np); %MPC Matrix Generation Function w/Input Penalization

%% Moving Window Blocking Setting
Nb=5; %Block Size
Nbi=1; %Blocking Index Counter

%% Simulation
Time=buoy.time(end)-Np*Ts; %Simulation Time
kT=round(Time/Ts); %Number of sample times

%Preallocation MPC Full DoF
X_FDoF=zeros(nx,kT); %MPC Full DoF States Matrix
Y_FDoF=zeros(ny,kT); %MPC Full DoF Outputs Matrix
U_FDoF=zeros(nu,kT); %MPC Full DoF Inputs Matrix
Energy_FDoF=zeros(1,kT); %MPC Full DoF Energy Absorbed Matrix
Power_FDoF=zeros(1,kT); %MPC Full DoF Generated Power Matrix
Timings_FDoF=zeros(1,kT); %MPC Full DoF Timings Matrix
Iterations_FDoF=zeros(1,kT); %MPC Full DoF Iterations Matrix

%Preallocation Standard Unconstrained MPC
X_FDoF_UNC=zeros(nx,kT); %Unconstrained MPC Full DoF States Matrix
Y_FDoF_UNC=zeros(ny,kT); %Unconstrained MPC Full DoF Outputs Matrix
U_FDoF_UNC=zeros(nu,kT); %Unconstrained MPC Full DoF Inputs Matrix
Energy_FDoF_UNC=zeros(1,kT); %Unconstrained MPC Full DoF Energy Absorbed Matrix
Power_FDoF_UNC=zeros(1,kT); %Unconstrained MPC Full DoF Generated Power Matrix

%Preallocation Standard UnconstrainedMPC + Penalisation on the control action
X_FDoF_UNC_LAM=zeros(nx,kT); %Unc. Penalized MPC Full DoF States Matrix
Y_FDoF_UNC_LAM=zeros(ny,kT); %Unc. Penalized MPC Full DoF Outputs Matrix
U_FDoF_UNC_LAM=zeros(nu,kT); %Unc. Penalized MPC Full DoF Inputs Matrix
Energy_FDoF_UNC_LAM=zeros(1,kT); %Unc. Penalized MPC Full DoF Energy Absorbed Matrix
Power_FDoF_UNC_LAM=zeros(1,kT); %Unc. Penalized MPC Full DoF Generated Power Matrix

%Preallocation MWB MPC
X_MWB=zeros(nx,kT); %MWB MPC States Matrix
Y_MWB=zeros(ny,kT); %MWB MPC Outputs Matrix
U_MWB=zeros(nu,kT); %MWB MPC Inputs Matrix
Energy_MWB=zeros(1,kT); %MWB MPC Energy Absorbed Matrix
Power_MWB=zeros(1,kT); %MWB MPC Generated Power Matrix
Timings_MWB=zeros(1,kT); %MWB MPC Timings Matrix
Iterations_MWB=zeros(1,kT); %MWB MPC Iterations Matrix

%Preallocation GPC
X_GPC=zeros(nx,kT); %GPC States Matrix
Y_GPC=zeros(ny,kT); %GPC Outputs Matrix
U_GPC=zeros(nu,kT); %GPC Inputs Matrix
Energy_GPC=zeros(1,kT); %GPC Energy Absorbed Matrix
Power_GPC=zeros(1,kT); %GPC Generated Power Matrix
Timings_GPC=zeros(1,kT); %GPC Timings Matrix
Iterations_GPC=zeros(1,kT); %GPC Iterations Matrix

%Quadprog Options Setup
options=optimoptions('quadprog','Display','off'); %Avoid Displaying Quadprog Messages

%Main For Loop
k_ini=3; %Starting Point of For Loop
for k=k_ini:kT-1
%Buoy Future Force
Future_Exc_Force=buoy.ExcForce(k:k+Np-1,1)/scaler; %Buoy Future Force Vector

%% Unconstrained MPC
%Compute augmented state
Xa_FDoF_UNC=[X_FDoF_UNC(:,k);
U_FDoF_UNC(:,k-1)/scaler;
(U_FDoF_UNC(:,k-1)-U_FDoF_UNC(:,k-2))/scaler]; %Augmented State

%Form linear term
f=f1*Xa_FDoF_UNC+f2*Future_Exc_Force; %Linear term
Uopt_FDoF_UNC=-H\f; %Unconstrained Solution

%Saturate Delta Solution
Uopt_FDoF_UNC(Uopt_FDoF_UNC<DeltaUmin)=DeltaUmin; %Saturate Lower DeltaU Limit
Uopt_FDoF_UNC(Uopt_FDoF_UNC>DeltaUmax)=DeltaUmax; %Saturate Upper DeltaU Limit

%Select first input
U_FDoF_UNC(:,k)=U_FDoF_UNC(:,k-1)+Uopt_FDoF_UNC(1:nu)*scaler; %Apply only first input

%Saturate Input
if U_FDoF_UNC(:,k)>Umax*scaler
U_FDoF_UNC(:,k)=Umax*scaler; %Saturate Upper U Limit
end
if U_FDoF_UNC(:,k)<Umin*scaler
U_FDoF_UNC(:,k)=Umin*scaler; %Saturate Lower U Limit
end

%System Simulation
X_FDoF_UNC(:,k+1) = Ad*X_FDoF_UNC(:,k) + Bd*U_FDoF_UNC(:,k) + Bd*buoy.ExcForce(k,1); %States
Y_FDoF_UNC(:,k) = Cd*X_FDoF_UNC(:,k); %Outputs
Energy_FDoF_UNC(:,k)=Energy_FDoF_UNC(:,k-1)-X_FDoF_UNC(2,k)*U_FDoF_UNC(1,k); %Energy Absorbed
Power_FDoF_UNC(:,k)=-X_FDoF_UNC(2,k)*U_FDoF_UNC(1,k); %Generated Power

%% Unconstrained MPC + Penalisation
%Compute augmented state
Xa_FDoF_UNC_LAM=[X_FDoF_UNC_LAM(:,k);
U_FDoF_UNC_LAM(:,k-1)/scaler;
(U_FDoF_UNC_LAM(:,k-1)-U_FDoF_UNC_LAM(:,k-2))/scaler]; %Augmented State

%Form linear term
f_lam=f1_lam*Xa_FDoF_UNC_LAM+f2_lam*Future_Exc_Force; %Linear term
Uopt_FDoF_UNC_LAM=-H_lam\f_lam; %Unconstrained Solution

%Saturate Optimal Solution
Uopt_FDoF_UNC_LAM(Uopt_FDoF_UNC_LAM<DeltaUmin)=DeltaUmin; %Saturate Lower DeltaU Limit
Uopt_FDoF_UNC_LAM(Uopt_FDoF_UNC_LAM>DeltaUmax)=DeltaUmax; %Saturate Upper DeltaU Limit

%Select first input
U_FDoF_UNC_LAM(:,k)=U_FDoF_UNC_LAM(:,k-1)+Uopt_FDoF_UNC_LAM(1:nu)*scaler; %Apply only first input

%Saturate Input
if U_FDoF_UNC_LAM(:,k)>Umax*scaler
U_FDoF_UNC_LAM(:,k)=Umax*scaler; %Saturate Upper U Limit
end
if U_FDoF_UNC_LAM(:,k)<Umin*scaler
U_FDoF_UNC_LAM(:,k)=Umin*scaler; %Saturate Lower U Limit
end

%System Simulation
X_FDoF_UNC_LAM(:,k+1) = Ad*X_FDoF_UNC_LAM(:,k) + Bd*U_FDoF_UNC_LAM(:,k) + Bd*buoy.ExcForce(k,1); %States
Y_FDoF_UNC_LAM(:,k) = Cd*X_FDoF_UNC_LAM(:,k); %Outputs
Energy_FDoF_UNC_LAM(:,k)=Energy_FDoF_UNC_LAM(:,k-1)-X_FDoF_UNC_LAM(2,k)*U_FDoF_UNC_LAM(1,k); %Energy Absorbed
Power_FDoF_UNC_LAM(:,k)=-X_FDoF_UNC_LAM(2,k)*U_FDoF_UNC_LAM(1,k); %Generated Power

%% Full MPC Optimization
%Compute augmented state
Xa_FDoF=[X_FDoF(:,k);
U_FDoF(:,k-1)/scaler;
(U_FDoF(:,k-1)-U_FDoF(:,k-2))/scaler]; %Augmented State

%Form Constraint matrices
M_inputdelta=[eye(Np);
-eye(Np);
tril(ones(Np));
-tril(ones(Np))]; %Constraints Matrix

%Standard quadprog
tic %Start Timer
f=f1*Xa_FDoF+f2*Future_Exc_Force; %Linear term
b_inputdelta=[repmat(DeltaUmax,Np,1);
-repmat(DeltaUmin,Np,1);
repmat(Umax,Np,1)-U_FDoF(:,k-1)/scaler;
-(repmat(Umin,Np,1)-U_FDoF(:,k-1)/scaler)]; %Input/Input Rate Constraints
Uopt_FDoF=quadprog(H,f,M_inputdelta,b_inputdelta,[],[],[],[],[],options); %Solve Quadprog
Timings_FDoF(:,k)=toc; %Save Time

%Capture iterations separately
[~,~,~,output]=quadprog(H,f,M_inputdelta,b_inputdelta,[],[],[],[],[],options); %Solve Quadprog for Iterations
Iterations_FDoF(:,k)=output.iterations; %Save Iterations

%Saturate Delta Optimal Solution
Uopt_FDoF(Uopt_FDoF<DeltaUmin)=DeltaUmin; %Saturate Lower DeltaU Limit
Uopt_FDoF(Uopt_FDoF>DeltaUmax)=DeltaUmax; %Saturate Upper DeltaU Limit

%Select first input
U_FDoF(:,k)=U_FDoF(:,k-1)+Uopt_FDoF(1:nu)*scaler; %Apply only first input

%Saturate Input
if U_FDoF(:,k)>Umax*scaler
U_FDoF(:,k)=Umax*scaler; %Saturate Upper U Limit
end
if U_FDoF(:,k)<Umin*scaler
U_FDoF(:,k)=Umin*scaler; %Saturate Lower U Limit
end

%System Simulation
X_FDoF(:,k+1) = Ad*X_FDoF(:,k) + Bd*U_FDoF(:,k) + Bd*buoy.ExcForce(k,1); %States
Y_FDoF(:,k) = Cd*X_FDoF(:,k); %Outputs
Energy_FDoF(:,k)=Energy_FDoF(:,k-1)-X_FDoF(2,k)*U_FDoF(1,k); %Energy absorbed
Power_FDoF(:,k)=-X_FDoF(2,k)*U_FDoF(1,k); %Generated Power

%% Moving Window Blocking MPC
%Calculate augmented state
Xa_MWB=[X_MWB(:,k);
U_MWB(:,k-1)/scaler;
(U_MWB(:,k-1)-U_MWB(:,k-2))/scaler]; %Augmented State

%Blocking Approach
N=Delta_Blocking(Np,Nb,Nbi,nu); %Moving Window Blocking Matrix
Hblk=N'*H*N; %Blocked/Compressed Hessian
Hblk=(Hblk+Hblk')/2; %Symmetrize Hessian
Nbt=size(N,2); %Number of decision variables
Mblk_inputdelta_red=[eye(Nbt);
-eye(Nbt);
tril(ones(Nbt));
-tril(ones(Nbt))]; %Constraint Matrix
f1_blk=N'*f1; %Blocked linear terms
f2_blk=N'*f2; %Blocked linear terms

%Measure times
tic %Start Timer
fblk=f1_blk*Xa_MWB+f2_blk*Future_Exc_Force; %Blocked linear term
b_inputdelta_red=[repmat(DeltaUmax,Nbt,1);
-repmat(DeltaUmin,Nbt,1);
repmat(Umax,Nbt,1)-U_MWB(:,k-1)/scaler;
-(repmat(Umin,Nbt,1)-U_MWB(:,k-1)/scaler)]; %Input/Input Rate Constraints
Uopt_MWB_blk=quadprog(Hblk,fblk,Mblk_inputdelta_red,b_inputdelta_red,[],[],[],[],[],options); %Solve Quadprog
Uopt_MWB=N*Uopt_MWB_blk; %De-compress
Timings_MWB(:,k)=toc; %Save Time

%Capture iterations separately
[~,~,~,output]=quadprog(Hblk,fblk,Mblk_inputdelta_red,b_inputdelta_red,[],[],[],[],[],options); %Solve Quadprog for Iterations
Iterations_MWB(:,k)=output.iterations; %Save Iterations

%Saturate Solution
Uopt_MWB(Uopt_MWB<DeltaUmin)=DeltaUmin; %Saturate Lower DeltaU Limit
Uopt_MWB(Uopt_MWB>DeltaUmax)=DeltaUmax; %Saturate Upper DeltaU Limit

%Select first input
U_MWB(:,k)=U_MWB(:,k-1)+Uopt_MWB(1:nu)*scaler; %Apply only first input

%Saturate Input
if U_MWB(:,k)>Umax*scaler
U_MWB(:,k)=Umax*scaler; %Saturate Upper U Limit
end
if U_MWB(:,k)<Umin*scaler
U_MWB(:,k)=Umin*scaler; %Saturate Lower U Limit
end

%System Simulation
X_MWB(:,k+1) = Ad*X_MWB(:,k) + Bd*U_MWB(:,k) + Bd*buoy.ExcForce(k,1); %States
Y_MWB(:,k) = Cd*X_MWB(:,k); %Outputs
Energy_MWB(:,k)=Energy_MWB(:,k-1)-X_MWB(2,k)*U_MWB(1,k); %Energy Absorbed
Power_MWB(:,k)=-X_MWB(2,k)*U_MWB(1,k); %Generated Power

%Increase block size and reset if at the limit
Nbi=Nbi+1; %Increase block size
if Nbi>Nb
Nbi=1; %Reset if block size reach limit
end

%% Generalized Predictive Control
%Compute augmented state
Xa_GPC=[X_GPC(:,k);
U_GPC(:,k-1)/scaler;
(U_GPC(:,k-1)-U_GPC(:,k-2))/scaler]; %Augmented State

%GPC Blocking Approach
N=Standard_Delta_Blocking(Np,Nb); %GPC Blocking Matrix
Hblk=N'*H*N; %Compressed Hessian
Hblk=(Hblk+Hblk')/2; %Symmetrize Hessian
Nbt=size(N,2); %Number of decision variables
Mblk_inputdelta_red=[eye(Nbt);
-eye(Nbt);
tril(ones(Nbt));
-tril(ones(Nbt))]; %Constraint Matrix
f1_blk=N'*f1; %Blocked linear terms
f2_blk=N'*f2;     %Blocked linear terms

%Measure times
tic
fblk=f1_blk*Xa_GPC+f2_blk*Future_Exc_Force; %Linear term
b_inputdelta_red=[repmat(DeltaUmax,Nbt,1);
-repmat(DeltaUmin,Nbt,1);
repmat(Umax,Nbt,1)-U_GPC(:,k-1)/scaler;
-(repmat(Umin,Nbt,1)-U_GPC(:,k-1)/scaler)]; %Input/Input Rate Constraints
Uopt_GPC_blk=quadprog(Hblk,fblk,Mblk_inputdelta_red,b_inputdelta_red,[],[],[],[],[],options); %Solve Quadprog
Uopt_GPC=N*Uopt_GPC_blk; %De-compress
Timings_GPC(:,k)=toc; %Save Time

%Capture iterations separately
[~,~,~,output]=quadprog(Hblk,fblk,Mblk_inputdelta_red,b_inputdelta_red,[],[],[],[],[],options); %Solve Quadprog for Iterations
Iterations_GPC(:,k)=output.iterations; %Save Iterations

%Saturate Solution
Uopt_GPC(Uopt_GPC<DeltaUmin)=DeltaUmin; %Saturate Lower DeltaU Limit
Uopt_GPC(Uopt_GPC>DeltaUmax)=DeltaUmax; %Saturate Upper DeltaU Limit

%Select first input
U_GPC(:,k)=U_GPC(:,k-1)+Uopt_GPC(1:nu)*scaler; %Apply only first input

%Saturate Input
if U_GPC(:,k)>Umax*scaler
U_GPC(:,k)=Umax*scaler; %Saturate Upper U Limit
end
if U_GPC(:,k)<Umin*scaler
U_GPC(:,k)=Umin*scaler; %Saturate Lower U Limit
end

%System Simulation
X_GPC(:,k+1) = Ad*X_GPC(:,k) + Bd*U_GPC(:,k) + Bd*buoy.ExcForce(k,1);     %States
Y_GPC(:,k) = Cd*X_GPC(:,k);    %Outputs
Energy_GPC(:,k)=Energy_GPC(:,k-1)-X_GPC(2,k)*U_GPC(1,k); %Energy Absorbed
Power_GPC(:,k)=-X_GPC(2,k)*U_GPC(1,k); %Generated Power

%% Plot Predictions
if k*Ts+1.1*Ts>=180&&k*Ts+1.1*Ts<180+Ts
Plot_All_Predictions(); %Plot Optimal Solution Predictions
end
end