并网逆变器控制策略——模型预测控制MPC（一）

下次一定毕业

已于 2024-06-19 17:23:41 修改

阅读量911

点赞数 31

分类专栏：三相并网逆变器 Simulink仿真文章标签： matlab 学习笔记经验分享算法

于 2024-05-08 21:12:44 首次发布

本文链接：https://blog.csdn.net/weixin_45933967/article/details/138545501

版权

三相并网逆变器同时被 2 个专栏收录

5 篇文章 0 订阅

订阅专栏

Simulink仿真

5 篇文章 15 订阅

订阅专栏

文章目录

前言
一、三相三线并网逆变器状态方程
二、构建代价函数
- 1.代价函数的定义形式
- 2.代价函数的二次规划形式
三、仿真效果
- 仿真结果
总结

前言

模型预测控制 (MPC) 是一种最优控制技术，在每个控制周期 tk , MPC 控制器获得系统（被控对象）当前的状态。接下来它利用基于系统当前状态和系统动态模型通过在有限时间上预测系统未来 p 步状态轨迹，并与目标轨迹构建代价函数和约束，进行一次开环优化问题（要优化的变量是作用在被控对象上的控制输入序列）求解，得到未来一段时间 [tk,tk+1,···,tk+p] 的控制输入序列 [uk,uk+1,···,uk+p] 。当然对于求解得到的控制输入序列，控制器只将最近时刻 tk 的控制 uk 作用于系统（被控对象）而忽略掉控制序列的其他值[uk+1,···,uk+p], 在下一个时刻tk+1，MPC 控制器会重复上述优化求解过程重新计算控制序输入列并只将第一个控制值作为当前时刻的控制量作用于系统，依次类推，重复上述过程进行下一时刻的求解。所以 MPC 在每个时刻都在线进行一个含约束的优化问题的求解（滚动优化，特殊情况除外）。
在这里插入图片描述
以并网逆变器为被控系统，需要构建逆变器的的状态方程以得到其代价函数，本文以三相三线制并网逆变器为对象，在不进行约束的情况下通过模型预测控制（MPC）进行控制。

一、三相三线并网逆变器状态方程

我们以三相三线制并网逆变器为控制对象，以三相对称的理想情况为例，如下图所示

在这里插入图片描述
其中 $e_{a}$ 、 $e_{b}$ 、 $e_{c}$ 为电网侧的电压，在理想状态下，电网稳定运行且容量无限大， $e_{a}$ 、 $e_{b}$ 、 $e_{c}$ 三相对称，幅值和频率均稳定不变。
通过基尔霍夫定律，我们可以得到在abc静止坐标系下的逆变器数学方程：

$L_{f} \begin{bmatrix} \frac{\mathrm{d}}{\mathrm{d}t}i_{a} \\ \frac{\mathrm{d}}{\mathrm{d}t}i_{b} \\ \frac{\mathrm{d}}{\mathrm{d}t}i_{c} \end{bmatrix} =\begin{bmatrix} U_{a} \\ U_{b} \\ U_{c} \end{bmatrix}-R_{f} \begin{bmatrix} i_{a}\\i_{b}\\ i_{c} \end{bmatrix}-\begin{bmatrix} e_{a}\\e_{b}\\ e_{c} \end{bmatrix}$

为简化控制，通过Park变换将其转换为dq旋转坐标系下的表达式：

$L_{f}\begin{bmatrix}\frac{\mathrm{d} i_{d}}{\mathrm{d} t} \\\frac{\mathrm{d} i_{q}}{\mathrm{d} t} \end{bmatrix} =\begin{bmatrix} U_{d} \\ U_{q} \end{bmatrix}-R_{f} \begin{bmatrix} i_{d}\\i_{q}\end{bmatrix}-\begin{bmatrix} e_{d}\\e_{q}\end{bmatrix}+\begin{bmatrix} \omega L_{f} i_{q}\\-\omega L_{f} i_{d}\end{bmatrix}$

其中 $U_{a}$ 、 $U_{b}$ 、 $U_{c}$ 为逆变器三相桥臂对应的电压， $i_{a}$ 、 $i_{b}$ 、 $i_{c}$ 为三相之路的电流， $U_{d}$ 、 $U_{q}$ 、 $i_{d}$ 、 $i_{q}$ 则分别表示对应dp坐标系下的分量。
到这里我们基本就得到了并网逆变器的数学模型，即状态方程表达式。但要通过MPC进行预测控制，需要对连续的状态方程进行离散化，以得到k时刻前后系统的输入输出关系式。这里我们采用欧拉离散的方法，即

$\dot{x}(k)=\frac{x(k+1)-x(k)}{T}$

其中 $T$ 为系统的采样时间，将该式代入dq坐标系下的状态方程并整理，即可得到离散后的系统状态方程：

$\begin{bmatrix}i_{d}(k+1) \\i_{q}(k+1) \end{bmatrix} =\begin{bmatrix}1-\frac{TR_{f}}{L_{f}} & wT\\ -wT &1-\frac{TR_{f}}{L_{f}}\end{bmatrix}\begin{bmatrix} i_{d}(k) \\ i_{q}(k) \end{bmatrix}+\begin{bmatrix} \frac{T}{L_{f}} &0 \\ 0 &\frac{T}{L_{f}}\end{bmatrix} \begin{bmatrix} U_{d}(k)-e_{d}(k)\\U_{q}(k)-e_{q}(k)\end{bmatrix}$

我们将该式整理为状态空间离散预测模型 $x (k + 1) = A x (k) + B u (k)$ 的形式，即 $x=\begin{bmatrix}i_{d} \\i_{q} \end{bmatrix}$ ； $u(k)=\begin{bmatrix} U_{d}(k)-e_{d}(k)\\U_{q}(k)-e_{q}(k)\end{bmatrix}$ ; $A=\begin{bmatrix}1-\frac{TR_{f}}{L_{f}} & wT\\ -wT &1-\frac{TR_{f}}{L_{f}}\end{bmatrix}$ ; $B=\begin{bmatrix} \frac{T}{L_{f}} &0 \\ 0 &\frac{T}{L_{f}}\end{bmatrix}$

二、构建代价函数

1.代价函数的定义形式

这一部分将针对并网逆变器的状态空间方程以及控制要求进行推导。系统的整体控制结构如下，其中 $U_{d}$ 、 $U_{q}$ 作为被控系统的输入，系统误差由 $i_{d}$ 、 $i_{q}$ 和 $i_{dref}$ 、 $i_{qref}$ 作差得到。
在这里插入图片描述
针对逆变器系统，代价函数由系统误差和输入的加权和构成，MPC控制器采取滚动优化计算过程，控制器通过求解代价函数的最小化问题，得到每轮预测区间的系统输入，以此过程循环滚动计算，最终使系统输出达到参考值。在上述过程中，我们得到了逆变器状态空间离散预测模型，以k时刻为起点，可以获得整个预测区间N步上的参考轨迹。
按代价函数的定义，系统误差为 $i_{d}-i_{dref}$ 、 $i_{q}-i_{qref}$ ，将其统称为 $i-i_{ref}$ 的向量,以 $E$ 来表示，通常代价函数有两种构成形式，一种以误差的绝对值形式构成，令一种则为平方形式，其目的都是为了使代价函数为大于等于0的数值，方便进行最小值计算，这里采用平方的形式，并以二次型的形式展现，则由系统误差和输入构成的代价函数为：

$J=\sum_{i=0}^{N-1} \left \{ \left [ i(k+i|k)-i_{ref} \right ] ^TQ\,\left [ i(k+i|k)-i_{ref} \right ]+u(k+i|k)^TR\,u(k+i|k)\right \}+\left [ i(k+N)-i_{ref} \right ] ^TF\,\left [ i(k+N)-i_{ref} \right ]$

$J=\sum_{i=0}^{N-1} \left \{ E(k+i|i) ^TQ\,E(k+i|i) +u(k+i|k)^TR\,u(k+i|k)\right \}+E(k+N) ^TF\,E(k+N)$

其中， $i (k + i ∣ i)$ 为在k时刻所预测的第i步的值， $N$ 为滚动计算的步长。 $\sum_{i=0}^{N-1} E(k+i|i) ^TQ\,E(k+i|i)$ 为系统的误差加权和， $\sum_{i=0}^{N-1} u(k+i|k)^TR\,u(k+i|k)$ 为系统输入加权和， $E(k+N) ^TF\,E(k+N)$ 为终端误差， $Q$ 、 $R$ 、 $F$ 分别为对应的权重系数，一般对对角矩阵。

2.代价函数的二次规划形式

代价函数的求解实际上可以转化为二次规划的问题，即 $min_{} J= X_{k}^{T}E\,u_{k}+u_{k}^{T}H\,u_{k}$ ,为了将上述代价函数转换为二次规划的形式，就需要利用到逆变器的状态空间离散预测模型，即:

$i (k + 1) = A i (k) + B u (k),$ $u (k) = U (k) - e (k)$

在此之前，需要对代价函数先进行进一步的整理。我们以 $N$ 为滚动计算的步数进行展开，则系统的每一步输入输出关系式可以表示为：

$\left\{\begin{matrix}E(k|k)=i_{k}-i_{ref} \\E(k+1|k)=i(k+1|k)-i_{ref}=Ai_{k}+Bu(k|k)-i_{ref} \\E(k+2|k)=i(k+2|k)-i_{ref}=Ai(k+1|k)+Bu(k+1|k)-i_{ref}=A^2i_{k}+ABu(k|k)-i_{ref} \\\vdots \\E(k+N|k)=A^Ni_{k}+A^{N-1}Bu(k|k)+\dots +Bu(k+N-1|k)-i_{ref}\end{matrix}\right.$

记 $E_{k}=\left\{\begin{matrix}E(k|k) \\E(k+1|k) \\E(k+2|k) \\\vdots \\E(k+N|k)\end{matrix}\right.$ ，则上述关系式可以以矩阵的形式表示：

$E_{k}=\begin{bmatrix}E \\A \\\vdots \\A^N\end{bmatrix}i_{k}+\begin{bmatrix} 0 & 0 & \cdots &0 \\ B& 0 & \cdots &0 \\ AB& B& \cdots & 0\\ \vdots & \vdots & &\vdots \\ A^{N-1}B& A^{N-2}B &\cdots &B\end{bmatrix}u_k-i_{ref}$

记 $M=\begin{bmatrix}E \\A \\\vdots \\A^N\end{bmatrix}$ , $C=\begin{bmatrix} 0 & 0 & \cdots &0 \\ B& 0 & \cdots &0 \\ AB& B& \cdots & 0\\ \vdots & \vdots & &\vdots \\ A^{N-1}B& A^{N-2}B &\cdots &B\end{bmatrix}$ , $u_k=\begin{bmatrix}u(k|k) \\u(k+1|k) \\\vdots \\u(k+N-1|k)\end{bmatrix}$ ，则表达式可以变为：

$E_k=Mi_k+Cu_k-i_{ref}$

其中 $i_k$ 为系统在k时刻的电流初值，为已知量， $u_k$ 为一次滚动计算的输入电压组成的向量。由上述推导过程，我们可以将代价函数做进一步变换；

$J=\sum_{i=0}^{N-1} \left \{ E(k+i|i) ^TQ\,E(k+i|i) +u(k+i|k)^TR\,u(k+i|k)\right \}+E(k+N) ^TF\,E(k+N)$

将代价函数中的误差加权和连同终端误差一起进行展开，可以写作 $\begin{bmatrix}E(k|k) \\E(k+1|k) \\\vdots \\E(k+N|k)\end{bmatrix}^T\begin{bmatrix} Q& & & \\ & Q & & \\ & &\ddots & \\ & & &F\end{bmatrix}\begin{bmatrix}E(k|k) \\E(k+1|k) \\\vdots \\E(k+N|k)\end{bmatrix}$ ；同理输入加权和也可以写作 $\begin{bmatrix}u(k|k) \\u(k+1|k) \\\vdots \\u(k+N-1|k)\end{bmatrix}^T\begin{bmatrix} R& & & \\ & R & & \\ & &\ddots & \\ & & &R\end{bmatrix}\begin{bmatrix}u(k|k) \\u(k+1|k) \\\vdots \\u(k+N-1|k)\end{bmatrix}$ 。由此可以将代价函数改写为：

$J=E_k^T\,\overline{Q} \,E_k+u_k^T\,\overline{R} \,u_k$

由于我们所要求解的是使代价函数最小的输入 $u_k$ ，而在这个表达式中还含有未知的 $E_k$ ，因此我们需要将 $E_k$ 以 $u_k$ 的形式来表达从而消去 $E_k$ ，将代价函数变为只含有 $u_k$ 的单变量函数，将 $E_k=Mi_k+Cu_k-i_{ref}$ 代入代价函数，即

$J=(Mi_k+Cu_k-i_{ref})^T\,\overline{Q} \,(Mi_k+Cu_k-i_{ref})+u_k^T\,\overline{R} \,u_k$

将其进一步展开，由于 $i_k$ 和 $i_{ref}$ 均为已知的固定值，因此在对代价函数进行最小求解时不会影响最终的解，所以只需要提取其中含有 $u_k$ 的项来对代价函数进行简化，得到：

$J=i_{k}^{T} M^T\,\overline{Q}\,C\,u_k+u_k^T\,C^T\,\overline{Q}\,M\,i_{k}+u_{k}^{T}\,C^T\,\overline{Q}\,C\,u_k-u_{k}^{T}\,C^T\,\overline{Q}\,i_{ref}-i_{ref}^T\,\overline{Q}\,C\,u_k+u_{k}^{T}\,\overline{R}u_{k}$

其中不难看出， $i_{k}^{T} M^T\,\overline{Q}\,C\,u_k$ 和 $u_k^T\,C^T\,\overline{Q}\,M\,i_{k}$ 互为转置的关系，而代价函数中的各项表达式均为一个一维数值，所以这两部分实际是相等的。同理， $u_{k}^{T}\,C^T\,\overline{Q}\,i_{ref}$ 和 $i_{ref}^T\,\overline{Q}\,C\,u_k$ 也是相等的，从而将代价函数进一步化简为：

$J=2i_{k}^{T} M^T\,\overline{Q}\,C\,u_k+u_{k}^{T}\,C^T\,\overline{Q}\,C\,u_k-2i_{ref}^T\,\overline{Q}\,C\,u_k+u_{k}^{T}\,\overline{R}u_{k}$

$J=(2i_{k}^{T} M^T\,\overline{Q}\,C-2i_{ref}^T\,\overline{Q}\,C\,)u_k+u_{k}^{T}\,(C^T\,\overline{Q}\,C+\overline{R}\,)\,u_k$

令 $Z=2i_{k}^{T} M^T\,\overline{Q}\,C-2i_{ref}^T\,\overline{Q}\,C$ ， $H=C^T\,\overline{Q}\,C+\overline{R}$ ，则代价函数最终可以化简为

$J=Zu_k+u_{k}^{T}H\,u_k$

该形式即可满足二次规划问题的求解形式，到这里，我们就将求解使代价函数最小的输入 $u_k$ 的问题转化为关于 $u_k$ 得二次规划问题，这样就可以通过二次规划的求解方式来得到 $u_k$ 使得代价函数最小，从而在一个滚动计算周期内得到合适的输入 $u_k$ ，使得系统误差达到最小，系统输入不断趋近于系统的参考值。

三、仿真效果

通过上述推导过程得到了代价函数的最终求解形式，以滚动计算步长N=5，迭代100次为例，通过Matlab进行验证。 $i_k$ 为系统的初值，由上述的逆变器状态方程， $i_k$ 由 $i_d$ 和 $i_q$ 构成，为1×2的向量，假定系统并网运行，输出电流为100A（幅值），没有无功注入，则 $i_{dref}=-100$ ， $i_{qref}=0$ 。以N=5的情况下， $i_{ref}$ 为12×1的向量，由 $i_{dref}$ 和 $i_{qref}$ 组成。设系统采样时间T=(1e-5)，Rf=0.05，L=0.002，默认各初值均为0

function [E , H, Z]=MPC_Matrices(A,B,Q,R,F,N)

n=size(A,1); % A 是 n x n 矩阵, 得到 n

p=size(B,2); % B 是 n x p 矩阵, 得到 p

M=[eye(n);zeros(N*n,n)]; % 初始化 M 矩阵. M 矩阵是 (N+1)n x n的， 

C=zeros((N+1)*n,N*p); % 初始化 C 矩阵, 这一步令它有 (N+1)n x NP 个 0

% 定义M 和 C 

tmp=eye(n); %定义一个n x n 的 I 矩阵

%　更新Ｍ和C

for i=1:N % 循环，i 从 1到 N

    rows =i*n+(1:n); %定义当前行数，从i x n开始，共n行 

    C(rows,:)=[tmp*B,C(rows-n, 1:end-p)]; %将c矩阵填满

    tmp= A*tmp; %每一次将tmp左乘一次A

    M(rows,:)=tmp; %将M矩阵写满

end
% 定义Q_bar和R_bar

Q_bar = kron(eye(N),Q);

Q_bar = blkdiag(Q_bar,F);

R_bar = kron(eye(N),R); 
% 计算G, E, H

E=M'*Q_bar*C; % E: NP x n

H=C'*Q_bar*C+R_bar; % NP x NP 

Z=Q_bar*C;

end

function u_k= Prediction(x_k,xref,E,H,N,Z,p)

U_k = zeros(N*p,1); % NP x 1
options = optimoptions('quadprog','Algorithm','active-set');
x0 = zeros(N*p,1);

U_k = quadprog((H+H')/2,2*x_k'*E-2*xref'*Z,[],[],[],[],[],[],x0,options);  % (H+H')/2 保证H矩阵对称

u_k = U_k(1:p,1); % 取第一个结果

end


%% 定义状态矩阵 A, n x n 矩阵

A = [1-(1e-5)*0.05/(0.002),100*pi*(1e-5);-100*pi*(1e-5),1-(1e-5)*0.05/(0.002)];

n= size (A,1);

%% 定义输入矩阵 B, n x p 矩阵

B = [(1e-5)/(0.002),0;0,(1e-5)/(0.002)];

p = size(B,2);

%% 定义Q矩阵，n x n 矩阵

Q=[400 0;0 700];

%% 定义F矩阵，n x n 矩阵

F=[1 0;0 1];

%% 定义R矩阵，p x p 矩阵

R=[1 0;0 1];

%% 定义step数量k

k_steps=100;
%% 定义预测区间K

N=5;

%% 定义矩阵 i_K， n x k 矩 阵

i_K = zeros(n,k_steps);
i_K(:,1) =[0;0];

%% 初始状态变量值， n x 1 向量
%iref = zeros((N+1)*2,1);
i_ref = -100;
iref = [i_ref;0;i_ref;0;i_ref;0;i_ref;0;i_ref;0;i_ref;0];

%% 定义输入矩阵 U_K， p x k 矩阵

%%V_K = zeros(p,k_steps);
%%e = [536.6;9.12];
U_K = zeros(p,k_steps);
%%U_K(:,1) = V_K(:,1) - e;

%% Call MPC_Matrices 函数 求得 E,H矩阵 

[E,H,Z]=MPC_Matrices(A,B,Q,R,F,N);

%% 计算每一步的状态变量的值

for k = 1 : k_steps

%% 求得U_K(:,k)

U_K(:,k) = Prediction(i_K(:,k),iref,E,H,N,Z,p);

%% 计算第k+1步时状态变量的值

i_K(:,k+1)=(A*i_K(:,k)+B*U_K(:,k));

end

仿真结果

在这里插入图片描述

总结

本文推导了模型预测控制在三相三线制并网逆变器的应用场景，并通过Matlab对推导过程进行了仿真，文中的代价函数推导过程及仿真代码来自于B站up主DR_Can，在其讲解基础之上对于并网逆变器的应用场景进行调整，最终得到本文内容，如有纰漏欢迎指正！
关于模型预测控制MPC的代价函数部分，由其定义所产生的形式并不是唯一的，而在并网逆变器这一应用场景下更是有着多种实现形式，本文仅给出了在dq静止坐标系下，以传统PI控制器的系统结构为基础进行的应用，更多内容会在后续的文章中进行讲解。文中使用的程序文件已上传，有需要自取。

下次一定毕业

关注

31
点赞
踩
28

收藏

觉得还不错? 一键收藏
打赏
0
评论
并网逆变器控制策略——模型预测控制MPC（一）

模型预测控制 (MPC) 是一种最优控制技术，在每个控制周期 tk , MPC 控制器获得系统（被控对象）当前的状态。接下来它利用基于系统当前状态和系统动态模型通过在有限时间上预测系统未来 p 步状态轨迹，并与目标轨迹构建代价函数和约束，进行一次开环优化问题（要优化的变量是作用在被控对象上的控制输入序列）求解，得到未来一段时间 [tk,tk+1,···,tk+p] 的控制输入序列 [uk,uk+1,···,uk+p] 。当然对于求解得到的控制输入序列，控制器只将最近时刻 tk 的控制 uk 作用于系统（
复制链接

扫一扫