CN104932262A - 基于mpc和pi控制方法的带电动汽车的微网调频方法 - Google Patents

Abstract

本发明采用了一种基于MPC和PI控制方法的带电动汽车的微网调频方法。本发明充分利用电动汽车车载电池本身作为一种储能装置的特性，并在电动汽车入网技术（V2G）的支持下，将其参与到微网频率调节中来以补偿微网孤网运行时自身负荷频率控制的不足。在本发明方案中，每一个控制周期，首先通过状态观测器得到当前时刻系统的状态量，采用MPC控制器进行求解最优控制输入变量及其相应的预测状态量，求解结果用来优化PI控制器中的参数，即重新设定比例参数和积分参数。最后，采用优化了控制参数的PI控制器实现对微网频率的调节控制，直达下一个控制周期的到来。本发明能够实现较好的控制效果、并具有全局鲁棒性能和较强的可靠性能。

Description

基于MPC和PI控制方法的带电动汽车的微网调频方法

技术领域

本发明涉及电力系统频率调节控制领域，尤其涉及一种基于MPC和PI控制方法的带电动汽车的微网调频控制策略

背景技术

微电网是一种将本地可再生能源发电系统、清洁能源发电系统、储能装置以及各类负载集成在一起的新型能源系统模式，正在受到越来越多的关注。微网系统可运行于并网和离网两种模式下，当微网运行于离网的孤立状态时必须要有自身的频率控制策略，以保证微网能够继续正常运行并为在重新并网作准备。当微网运行在孤网状态时，功率不平衡时的频率控制以及不平衡功率的分配问题成为研究焦点，因为经过逆变器接入微网的微源(比如说风电、光伏)对频率变化不会自动产生响应，而像小水电、微型汽轮机等响应速度比较慢，而且具有较小的惯性常数。传统电网中当有冲击负荷出现时，可以由系统的惯性响应实现减少频率波动的目的，而微网不能依靠或者不能完全依靠发电机的惯性来实现频率波动的减少，电动汽车车载电池本身作为一种储能装置，在电动汽车入网技术(V2G)的支持下，已经实现了与电网之间电能的双向互动。随着电动汽车入网技术以及微网技术的成熟，将电动汽车纳入微网系统，充分利用车载电池存储的电能来参与微网调频以补偿微网自身负荷频率控制的不足。因此，一个有效的协调控制策略是实现V2G参与微网调频的关键。

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种没有改变原有PI控制系统的结构，即可使原有基于PI控制的高级策略仍能正常发挥作用，具有更好的兼容性的基于MPC和PI控制方法的带电动汽车的微网调频方法。

本发明还有一目的是解决现有技术所存在的技术问题；提供了一种可以最大限度在一个控制周期内让_PI控制器接近最优的基于MPC和PI控制方法的带电动汽车的微网调频方法。

本发明再有一目的是解决现有技术所存在的技术问题；提供了一种克服了PI控制器对运行工况的适应性，动态性能也得到积极改善；可以有效地避免MPC问题求解不可行时无法进行控制的问题，这时整个系统仍可有效运行，加强了控制方法的可靠性的基于MPC和PI控制方法的带电动汽车的微网调频方法。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种基于MPC和PI控制方法的带电动汽车的微网调频方法，其特征在于，包括以下步骤：

步骤1：检测被控微网的各个组件的系统参数，包括各电动汽车充放电站的系统参数并分析各个微网组件的负荷频率控制特性，确定系统状态量，并根据检测到的系统参数得到考虑线性连续时间系统的状态空间模型，该状态空间模型基于以下公式：

$\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=A_{c}x\left(t\right)+B_{cu}u\left(t\right)+B_{cd}d\left(t\right)\\ y_c\left(t\right)=C_{c}x\left(t\right)\end{array}\right.$     式一

将连续系统转化为线性离散时间系统，则线性离散时间系统的状态空间模型可表示如下：

$\left\{\begin{array}{c}x(k+1)=Ax\left(k\right)+B_{u}u\left(k\right)+B_{d}d\left(k\right)\\ y_c\left(k\right)=C_{c}x\left(k\right)\end{array}\right.$     式二

其中是状态变量；是控制输入变量；是被控输出变量；

其中 $A=e^{A_c{T}_s},B_u={\∫}_0^{Ts}{e}^{A_c\mathrm{\τ}}d\mathrm{\τ}\·B_{cu},B_d={\∫}_0^{T_s}{e}^{A_c\mathrm{\τ}}d\mathrm{\τ}\·B_d,$ T_s是系统采样时间；

为了引入积分以减少或消除静态误差，将式二改为增量模型：

$\left\{\begin{array}{c}\mathrm{\Δ}x(k+1)=A\mathrm{\Δ}x\left(k\right)+B_u\mathrm{\Δ}u\left(k\right)+B_d\mathrm{\Δ}d\left(k\right)\\ y_c\left(k\right)=C_c\mathrm{\Δ}x\left(k\right)+y_c(k-1)\end{array}\right.$     式三

其中， $\begin{array}{c}\mathrm{\Δ}x\left(k\right)=x\left(k\right)-x(k-1),\\ \mathrm{\Δ}u\left(k\right)=u\left(k\right)-u(k-1),\\ \mathrm{\Δ}x\left(k\right)=d\left(k\right)-d(k-1).\end{array}$

步骤2：采用MPC控制器读取当前控制周期内的状态量，进行预测控制，求解预测时域内最优控制变量，具体包括：

步骤2.1：状态监测，PMU测量系统的实时状态信息，即各状态量的值，由MGDS进行收集整理，并将数据发送至MPC控制器；

步骤2.2：预测未来；根据MGDS收集处理的实时状态信息，MPC控制器根据状态空间方程预测输出偏差Δy的动态轨迹；并根据建立控制量和预测量之间的关系，为最优控制决策奠定基础；

步骤2.3：最优控制决策，基于步骤2的预测输出结果，MPC控制器根据目标函数以及约束条件求解最优控制量；

系统输出的参考轨迹一般是可以确定的；由于扰动等因素，系统实际输出和预测结果有误差，可定义二次型性能指标函数，找出满足J最小的最优控制输入序列Δu(k|k),…,Δu(k+m-1|k)；

$\underset{\mathrm{\Δ}u\left(k\right|k),...,\mathrm{\Δ}u\left(k+m-1|k\right)}{\min }J=min\underset{i=1}{\overset{p}{\Σ}}{q}_i{\left(\hat{y}\right(k+i\left|k\right)-y_r\left(k+i|k\right))}^2+\underset{j=1}{\overset{n}{\Σ}}{r}_j{\left(\mathrm{\Δ}u\right(k+j-1\left|k\right))}^2$    式十六

其中，q_i和r_i是权重系数；yr(k+i|k)为系统输出参考轨迹；若将用向量表示为ΔU(k|k).

则可将式十一写成：

J＝[E(k)+S_uΔU(k|k)]^TQ[E(k)+S_uΔU(k|k)]+ΔU(k|k)^TRΔU(k|k)  式十七

其中，E(k)＝S_xΔx(k)+Γy_c(k)+S_dΔd(k)，

然后将式十七展开得到：

J＝E(k)^TQE(k)+E(k)^TQSΔU(k|k)+ΔU(k|k)^TS^TQE(k)+ΔU(k|k)^TS^TQSΔU(k|k)+ΔU(k|k)^TRΔU(k|k)   式十八

注意到E(k)^TQSΔU(k|k)＝ΔU(k|k)^TS^TQE(k)，而且E(k)^TQE(k)这一项与控制输入量无关，则式十八可写成：

$\underset{\mathrm{\Δ}U\left(k\right|k)}{\min }J=\underset{\mathrm{\Δ}U\left(k\right|k)}{\min }\frac{1}{2}\mathrm{\Δ}U{\left(k\right|k)}^{T}H\mathrm{\Δ}U\left(k\right|k)+f{\left(k\right)}^T\mathrm{\Δ}U\left(k\right|k)$    式十九其中，H＝2[S^TQS+R]，f(k)＝2S^TQE(k)，可把系统状态量和系统输出量都看成是系统控制输入的函数，由系统控制输入决定，定义U＝[u(k)u(k+1)…u(k+N)]^T，则上述问题也就变为：

$\begin{array}{c}J\left(k\right)=min\frac{1}{2}{U}^{T}HU+h^{T}U\\ \begin{array}{cc}s.t.& GU\≤g\end{array}\end{array}$       式二十

其中，G为不等式约束矩阵线性，最优控制决策就是对上式二十的二次规划问题的求解；

步骤3：根据步骤2求解到的最优控制变量及其相应的预测状态量，采用最小二乘法求解当前控制时域内的PI控制器的最佳控制系数，比例系数和积分系数并更新PI控制器的控制系数；

考虑到PI控制器有如下控制规律：

$u=K_p\stackrel{\·}{\mathrm{\ρ}}+K_i\mathrm{\ρ}$        式四

$\mathrm{\ρ}=\underset{0}{\overset{\mathrm{\∞}}{\∫}}(y-r)dt=\underset{0}{\overset{\mathrm{\∞}}{\∫}}(Cx-r)dt$        式五

其中，r为系统输出参考量；显然有：

$\stackrel{\·}{\mathrm{\ρ}}=Cx-r$         式六

则可根据步骤2得到的最优解 $U_t^{*}={\left[\begin{array}{cccc}{U}_{t_0}& U_{t_1}& ...& U_{t_{N-1}}\end{array}\right]}^T,$ 求解其对应的状态量 ${\left[\begin{array}{c}x\\ \mathrm{\ρ}\end{array}\right]}_t^{*}=\left[\begin{array}{cccc}{x}_{t_0}& x_{t_1}& ...& x_{t_{N-1}}\\ {\mathrm{\ρ}}_{t_0}& {\mathrm{\ρ}}_{t_1}& ...& {\mathrm{\ρ}}_{t_{N-1}}\end{array}\right],$ 通过最小二乘法求取当前预测时域内_T内最佳的比例、积分系数；

步骤4：在当前控制周期内，采用更新了控制系数的PI控制器实现对微网频率的调节控制。

在上述的一种基于MPC和PI控制方法的带电动汽车的微网调频方法，所述步骤2.2的具体方法是：根据预测控制的基本原理，首先以步骤2.1得到的最新测量值为初始条件，基于式三预测系统未来的动态；为此，设定预测时域为p，控制时域为m且m≤p；并定义：

定义一：控制时域之外，控制变量不变，即Δu(k+i)＝0,i＝m,m+1,…,p-1；

定义二：可测干扰在k时刻之后不变，即Δd(k+i)＝0,i＝1,2,…,p-1；

定义一是因为控制时域有可能小于预测时域，而预测系统未来动态需要在整个预测时域的控制输入；定义二是因为在当前k时刻，尚未获取干扰的未来取值；

在当前时刻k，测量值为x(k)，计算Δx(k)＝x(k)-x(k-1)；这个Δx(k)将作为预测系统未来动态的起点；由式三预测k+1到k+3时刻的状态增量如下：

Δx(k+1|k)＝AΔx(k)+B_uΔu(k)+B_dΔd(k)              式七

Δx(k+2|k)＝AΔx(k+1|k)+B_uΔu(k+1)+B_dΔd(k+1)

                                          式八

＝A²Δx(k)+AB_uΔu(k)+B_uΔu(k+1)+AB_dΔd(k)

Δx(k+3|k)＝AΔx(k+2|k)+B_uΔu(k+2)+B_dΔd(k+2)

                                             式九

＝A³Δx(k)+A²B_uΔu(k)+AB_uΔu(k+1)+B_uΔu(k+2)+A²B_dΔd(k)

上式中，k+1|k表示k时刻对k+1时刻的预测；符号“|”后面的k表示当前时刻为k；进而，可以预测k+m至k+p时刻的状态

  式十

进一步，有输出方程(2)可以预测k+1至k+p的被控输出

y_c(k+1|k)＝C_cΔx(k+1|k)+y_c(k)＝C_cAΔx(k)+C_cB_uΔu(k)+C_cB_dΔd(k)+y_c(k)  式十一

y_c(k+2|k)＝C_cΔx(k+2|k)+y_c(k+1|k)

＝(C_cA²+C_cA)Δx(k)+(C_cAB_u+C_cB_u)Δu(k)+C_cB_uΔu(k+1)+(C_cAB_d+C_cB_d)Δd(k)+y_c(k)  式十二

  式十三

$\begin{array}{c}{y}_c(k+p|k)=C_c\mathrm{\Δ}x(k+p\left|k\right)+y_c(k+p-1|k)=\underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^i\mathrm{\Δ}x\left(k\right)+\underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^{i-1}{B}_u\mathrm{\Δ}u\left(k\right)+\underset{i=1}{\overset{p-1}{\Σ}}{C}_c{A}^{i-1}{B}_u\mathrm{\Δ}u(k+1)\\ +...+\underset{i=1}{\overset{p-m+1}{\Σ}}{C}_c{A}^{i-1}{B}_u\mathrm{\Δ}u\left(k+m-1\right)+\underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^{i-1}{B}_d\mathrm{\Δ}d\left(k\right)+y_c\left(k\right),\end{array}$ 式十四定义p步预测输出向量和m输入向量如下：

$Y_p(k+1|k)\stackrel{def}{=}{\left[\begin{array}{c}{y}_c(k+1|k)\\ y_c(k+2|k)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_c(k+p|k)\end{array}\right]}_{p\×1},\mathrm{\Δ}U\left(k\right)\stackrel{def}{=}{\left[\begin{array}{c}\mathrm{\Δ}u\\ \mathrm{\Δ}u\left(k+1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \mathrm{\Δ}u\left(k+m-1\right)\end{array}\right]}_{p\×1}$

那么，对系统未来p步预测的输出可以由下面的预测方程计算：

Y_p(k+1|k)＝S_xΔx(k)+Γy_c(k)+S_dΔd(k)+S_uΔU(k)      式十五

其中：

因此，本发明具有如下优点：1.采用的这种方法同时结合了MPC控制方式和PI控制方式，这种方法比单纯的MPC和PI控制方式有更多的优势。2.它没有改变原有PI控制系统的结构，即可使原有基于PI控制的高级策略仍能正常发挥作用，即意味着更好的兼容性；3.可以最大限度在一个控制周期内让PI控制器接近最优，4.控制周期是一直向前推进的，因此PI控制器的参数也是在不断更新，克服了PI控制器对运行工况的适应性，动态性能也得到积极改善；5同时由于采用了有关判断逻辑，可以有效地避免MPC问题求解不可行时无法进行控制的问题，这时整个系统仍可有效运行，加强了控制方法的可靠性。

附图说明

图1是本发明涉及的MPC与PI整合示意图。

图2是本发明涉及的MPC-PC的流程图。

图3是本发明涉及的微网系统。

图4是本发明涉及的微网频率控制体系。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

本发明在如图3所示的微网系统体系中实施了基于GPC和PI控制方法的带电动汽车的微网调频方法。图3中微网系统的状态空间模型如下式所示：

$\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=A_{c}x\left(t\right)+B_{cu}u\left(t\right)+B_{cd}d\left(t\right)\\ y_c\left(t\right)=C_{c}x\left(t\right)\end{array}\right.$        (1)

其中，x＝[Δf ΔP_MT ΔX_MT ΔP_E1 ΔP_E2]^T，u＝[Δu_MT Δu_E1 Δu_E2]^T，w＝[ΔP_D]^T，y＝[Δf]^T。

$A=\left[\begin{array}{ccccc}0& \frac{1}{H_t}& 0& \frac{1}{H_t}& \frac{1}{H_t}\\ 0& -\frac{1}{T_t}& \frac{1}{T_t}& 0& 0\\ -\frac{1}{{\mathrm{RT}}_f}& 0& -\frac{1}{T_f}& 0& 0\\ 0& 0& 0& -\frac{1}{T_{e1}}& 0\\ 0& 0& 0& 0& -\frac{1}{T_{e2}}\end{array}\right],B={\left[\begin{array}{ccccc}0& 0& \frac{1}{T_f}& 0& 0\\ 0& 0& 0& \frac{1}{T_{e1}}& 0\\ 0& 0& 0& 0& \frac{1}{T_{e2}}\end{array}\right]}^T$

$D={\left[\begin{array}{ccccc}-\frac{1}{H_t}& 0& 0& 0& 0\end{array}\right]}^T,$ C＝[1 0 0 0 0]。

A，B，C，D分别是系统的状态矩阵，输入矩阵，扰动矩阵和输出矩阵。x，u，w分别是系统的状态变量、可控制输入变量和不可控输入变量。T_e1是电动汽车充放电站1的时间常数，Δu_E1是电动汽车充放电站1接收的负荷频率控制信号；T_e2是电动汽车充放电站2的时间常数，Δu_E2是电动汽车充放电站2接收的负荷频率控制信号；Δf是频率偏差，Δu_MT是微型燃气轮机接收的负荷频率控制信号，ΔX_MT是燃料系统阀门位置变化的增量，T_f是燃料系统的时间常数，T_t是燃气涡轮的时间常数，R是调速器系数，ΔP_D是扰动功率且由负荷扰动功率ΔP_L和风电系统扰动功率ΔP_w组成，H_t代表微网的惯性常数。

以下结合附图对本发明作进一步的详细描述。在每一个控制周期都执行以下步骤，具体的控制流程如图2所示。

步骤1：状态监测，获取测量单元传送过来的系统输出测量值，计算出与给定控制目标的偏差e，对e的大小进行判断，如果当前偏差小于设定的限制，进入步骤2；否则进入步骤3；

步骤2：当前的PI控制器的控制参数不作修改，继续采用该PI控制器实现对微网频率的调节控制，直到该控制周期结束，到下一个控制周期，回到步骤1.

步骤3：MGDC进行收集和整理由PMU测量系统测量到的实时状态信息，包括Δf，ΔP_E1，ΔP_E2，ΔP_MT和ΔP_D，并将数据发送至MPC控制器。

步骤4：预测未来。根据MGDS收集处理的实时状态信息，MPC控制器根据状态空间方程预测输出偏差Δy的动态轨迹。并根据建立控制量和预测量之间的关系，为最优控制决策奠定基础。

步骤4.1：设定预测时域为p，控制时域为m且m≤p，并作如下的假设：假设一：控制时域之外，控制变量不变，即Δu(k+i)＝0,i＝m,m+1,…,p-1

假设二：可测干扰在k时刻之后不变，即Δd(k+i)＝0,i＝1,2,…,p-1

步骤4.2根据步骤3得到的当前时刻k的状态量测量值，x(k)，计算Δx(k)＝x(k)-x(k-1)，将这个Δx(k)将作为预测系统未来动态的起点。预测k+1到k+p时刻的状态(实际上是状态增量，简单起见称为状态)如下：

Δx(k+1|k)＝AΔx(k)+B_uΔu(k)+B_dd(k),      (3)

Δx(k+2|k)＝AΔx(k+1|k)+B_uΔu(k+1)+B_dΔd(k+1)

                                         (4)

＝A²Δx(k)+AB_uΔu(k)+B_uΔu(k+1)+AB_dΔd(k)

  (5)

Δx(k+m|k)＝AΔx(k+m-1|k)+B_uΔu(k+m-1)+B_dΔd(k+m-1)

                                                      (6)

＝A^mΔx(k)+A^m-1B_uΔu(k)+A^m-2B_uΔu(k+1)+…+B_uΔu(k+m-1)+A^m-1B_dΔd(k),

Δx(k+p|k)＝AΔx(k+p-1|k)+B_uΔu(k+p-1)+B_dΔd(k+p-1)

                                                       (7)

＝A^pΔx(k)+A^p-1B_uΔu(k)+A^p-2B_uΔu(k+1)+…+A^p-mB_uΔu(k+m-1)+A^p-1B_dΔd(k),

上式中，k+1|k表示k时刻对k+1时刻的预测；符号“|”后面的k表示当前时刻为k。

步骤4.3：预测k+1至k+p的被控输出量。

y_c(k+1|k)＝C_cΔx(k+1|k)+y_c(k)＝C_cAΔx(k)+C_cB_uΔu(k)+C_cB_dΔd(k)+y_c(k)  (8)

y_c(k+2|k)＝(C_cA²+C_cA)Δx(k)+(C_cAB_u+C_cB_u)Δu(k)+C_cB_uΔu(k+1)+(C_cAB_d+C_cB_d)Δd(k)+y_c(k)  (10)

  (11)

  (12)

定义p步预测输出向量和m输入向量如下：

$Y_p(k+1|k)\stackrel{def}{=}{\left[\begin{array}{c}{y}_c(k+1|k)\\ y_c(k+2|k)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_c(k+p|k)\end{array}\right]}_{p\×1},\mathrm{\Δ}U\left(k\right)\stackrel{def}{=}{\left[\begin{array}{c}\mathrm{\Δ}u\\ \mathrm{\Δ}u\left(k+1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \mathrm{\Δ}u\left(k+m-1\right)\end{array}\right]}_{p\×1}$

那么，对系统未来p步预测的输出可以由下面的预测方程计算：

Y_p(k+1|k)＝S_xΔx(k)+Γy_c(k)+S_dΔd(k)+S_uΔU(k)   (13)

其中：

步骤5：最优控制决策，基于步骤4，协调控制器根据目标函数以及约束条件求解最优控制量。如若目标函数有解，进入步骤6；否则进入步骤7

系统输出的参考轨迹一般是可以确定的；由于扰动等因素，系统实际输出和预测结果有误差，可定义二次型性能指标函数，找出满足J最小的最优控制输入序列Δu(k|k),…,Δu(k+m-1|k)。

$\underset{\mathrm{\Δ}u\left(k\right|k),...,\mathrm{\Δ}u\left(k+m-1|k\right)}{\min }J=min\underset{i=1}{\overset{p}{\Σ}}{q}_i{\left(\hat{y}\right(k+i\left|k\right)-y_r\left(k+i|k\right))}^2+\underset{j=1}{\overset{n}{\Σ}}{r}_j{\left(\mathrm{\Δ}u\right(k+j-1\left|k\right))}^2---\left(14\right)$

其中，q_i和r_i是权重系数；y_r(k+i|k)为系统输出参考轨迹；若将Δu(k|k),…,Δu(k+m-1|k)用向量表示为ΔU(k|k).

则可将式(11)写成：

J＝[E(k)+S_uΔU(k|k)]^TQ[E(k)+S_uΔU(k|k)]+ΔU(k|k)^TRΔU(k|k)        (15)

其中，E(k)＝S_xΔx(k)+Γy_c(k)+S_dΔd(k)，

然后将式(15)展开得到：

J＝E(k)^TQE(k)+E(k)^TQSΔU(k|k)+ΔU(k|k)^TS^TQE(k)+ΔU(k|k)^TS^TQSΔU(k|k)+ΔU(k|k)^TRΔU(k|k)  (16)

注意到E(k)^TQSΔU(k|k)＝ΔU(k|k)^TS^TQE(k)，而且E(k)^TQE(k)这一项与控制输入量无关，则式(16)可写成：

$\underset{\mathrm{\Δ}U\left(k\right|k)}{\min }J=\underset{\mathrm{\Δ}U\left(k\right|k)}{\min }\frac{1}{2}\mathrm{\Δ}U{\left(k\right|k)}^{T}H\mathrm{\Δ}U\left(k\right|k)+f{\left(k\right)}^T\mathrm{\Δ}U\left(k\right|k)---\left(17\right)$

其中，H＝2[S^TQS+R]，f(k)＝2S^TQE(k)

可把系统状态量和系统输出量都看成是系统控制输入的函数，由系统控制输入决定，定义U＝[u(k)u(k+1)…u(k+N)]^T，则上述问题也就变位：

$\begin{array}{c}J\left(k\right)=min\frac{1}{2}{U}^{T}HU+h^{T}U\\ \begin{array}{cc}s.t.& GU\≤g\end{array}\end{array}---\left(18\right)$

其中，G为不等式约束矩阵线性.

最优控制决策就是对上式(18)的二次规划问题的求解。

步骤6：修正PI控制器的控制参数，比例参数K_p和积分参数K_i。利用最小二乘法求出在当前控制内T内的最佳比例、积分系数；然后对PI控制器的比例、积分系数进行更新。采用更新了控制参数的PI控制器实现对微网的频率进行调节，直到当前控制周期结束，到下一控制周期，回到步骤1.

考虑到PI控制器有如下控制规律：

$u=K_p\stackrel{\·}{\mathrm{\ρ}}+K_i\mathrm{\ρ}---\left(19\right)$

$\mathrm{\ρ}=\underset{0}{\overset{\mathrm{\∞}}{\∫}}(y-r)dt=\underset{0}{\overset{\mathrm{\∞}}{\∫}}(Cx-r)dt---\left(20\right)$

其中，r为系统输出参考量。显然有：

$\stackrel{\·}{\mathrm{\ρ}}=Cx-r---\left(21\right)$

则可根据步骤5得到的最优解 $U_t^{*}={\left[\begin{array}{cccc}{U}_{t_0}& U_{t_1}& ...& U_{t_{N-1}}\end{array}\right]}^T$ 并求出其对应的状态量 ${\left[\begin{array}{c}x\\ \mathrm{\ρ}\end{array}\right]}_t^{*}=\left[\begin{array}{cccc}{x}_{t_0}& x_{t_1}& ...& x_{t_{N-1}}\\ {\mathrm{\ρ}}_{t_0}& {\mathrm{\ρ}}_{t_1}& ...& {\mathrm{\ρ}}_{t_{N-1}}\end{array}\right],$ 通过最小二乘法求取当前控制时域内T内最佳的比例、积分系数。

步骤7：当前的PI控制器的控制参数不作修改，继续采用该PI控制器实现对微网频率的调节控制，直到该控制周期结束，到下一个控制周期，回到步骤1.

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (2)

1.一种基于MPC和PI控制方法的带电动汽车的微网调频方法，其特征在于，包括以下步骤：

步骤1：检测被控微网的各个组件的系统参数，包括各电动汽车充放电站的系统参数并分析各个微网组件的负荷频率控制特性，确定系统状态量，并根据检测到的系统参数得到考虑线性连续时间系统的状态空间模型，该状态空间模型基于以下公式：

$\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=A_{c}x\left(t\right)+B_{cu}u\left(t\right)+B_{cd}d\left(t\right)\\ y_c\left(t\right)=C_{c}x\left(t\right)\end{array}\right.$   式一

将连续系统转化为线性离散时间系统，则线性离散时间系统的状态空间模型可表示如下：

$\left\{\begin{array}{c}x(k+1)=Ax\left(k\right)+B_{u}u\left(k\right)+B_{d}d\left(k\right)\\ y_c\left(k\right)=C_{c}x\left(k\right)\end{array}\right.$   式二

其中是状态变量；是控制输入变量；是被控输出变量；

其中 $A=e^{A_c{T}_s},B_u={\∫}_0^{Ts}{e}^{A_c\mathrm{\τ}}d\mathrm{\τ}\·B_{cu},B_d={\∫}_0^{T_s}{e}^{A_c\mathrm{\τ}}d\mathrm{\τ}\·B_d,$ T_s是系统采样时间；

为了引入积分以减少或消除静态误差，将式二改为增量模型：

$\left\{\begin{array}{c}\mathrm{\Δ}x(k+1)=A\mathrm{\Δ}x\left(k\right)+B_u\mathrm{\Δ}u\left(k\right)+B_d\mathrm{\Δ}d\left(k\right)\\ y_c\left(k\right)=C_c\mathrm{\Δ}x\left(k\right)+y_c(k-1)\end{array}\right.$   式三

其中， $\begin{array}{c}\mathrm{\Δ}x\left(k\right)=x\left(k\right)-x(k-1),\\ \mathrm{\Δ}u\left(k\right)=u\left(k\right)-u(k-1),\\ \mathrm{\Δ}x\left(k\right)=d\left(k\right)-d(k-1).\end{array}$

步骤2：采用MPC控制器读取当前控制周期内的状态量，进行预测控制，求解预测时域内最优控制变量，具体包括：

步骤2.1：状态监测，PMU测量系统的实时状态信息，即各状态量的值，由MGDS进行收集整理，并将数据发送至MPC控制器；

步骤2.2：预测未来；根据MGDS收集处理的实时状态信息，MPC控制器根据状态空间方程预测输出偏差Δy的动态轨迹；并根据建立控制量和预测量之间的关系，为最优控制决策奠定基础；

步骤2.3：最优控制决策，基于步骤2的预测输出结果，MPC控制器根据目标函数以及约束条件求解最优控制量；

系统输出的参考轨迹一般是可以确定的；由于扰动等因素，系统实际输出和预测结果有误差，可定义二次型性能指标函数，找出满足J最小的最优控制输入序列Δu(k|k),…,Δu(k+m-1|k)；

$\underset{\mathrm{\Δ}u\left(k\right|k),\mathrm{...},\mathrm{\Δ}u(k+m-1|k)}{\min }J=min\underset{i=1}{\overset{p}{\Σ}}{q}_i{\left(\hat{y}\right(k+i\left|k\right)-y_r(k+i|k\left)\right)}^2+\underset{j=1}{\overset{n}{\Σ}}{r}_j{\left(\mathrm{\Δ}u\right(k+j-1\left|k\right))}^2$   式十六

其中，q_i和r_i是权重系数；y_r(k+i|k)为系统输出参考轨迹；若将Δu(k|k),…,Δu(k+m-1|k)用向量表示为ΔU(k|k).

则可将式十一写成：

J＝[E(k)+S_uΔU(k|k)]^TQ[E(k)+S_uΔU(k|k)]+ΔU(k|k)^TRΔU(k|k)  式十七

其中，E(k)＝S_xΔx(k)+Γy_c(k)+S_dΔd(k)，

然后将式十七展开得到：

J＝E(k)^TQE(k)+E(k)^TQSΔU(k|k)+ΔU(k|k)^TS^TQE(k)+ΔU(k|k)^TS^TQSΔU(k|k)+ΔU(k|k)^TRΔU(k|k)  式十八

注意到E(k)^TQSΔU(k|k)＝ΔU(k|k)^TS^TQE(k)，而且E(k)^TQE(k)这一项与控制输入量无关，则式十八可写成：

$\underset{\mathrm{\Δ}U\left(k\right|k)}{\min }J=\underset{\mathrm{\Δ}U\left(k\right|k)}{min}\frac{1}{2}\mathrm{\Δ}U{\left(k\right|k)}^{T}H\mathrm{\Δ}U\left(k\right|k)+f{\left(k\right)}^T\mathrm{\Δ}U\left(k\right|k)$   式十九

其中，H＝2[S^TQS+R]，f(k)＝2S^TQE(k)，可把系统状态量和系统输出量都看成是系统控制输入的函数，由系统控制输入决定，定义U＝[u(k) u(k+1) … u(k+N)]^T，则上述问题也就变为：

$J\left(k\right)=min\frac{1}{2}{U}^{T}HU+h^{T}U$   式二十

s.t.GU≤g

其中，G为不等式约束矩阵线性，最优控制决策就是对上式二十的二次规划问题的求解；

步骤3：根据步骤2求解到的最优控制变量及其相应的预测状态量，采用最小二乘法求解当前控制时域内的PI控制器的最佳控制系数，比例系数和积分系数并更新PI控制器的控制系数；

考虑到PI控制器有如下控制规律：

$u=K_p\stackrel{\·}{\mathrm{\ρ}}+K_i\mathrm{\ρ}$   式四

$\mathrm{\ρ}=\underset{0}{\overset{\mathrm{\∞}}{\∫}}(y-r)dt=\underset{0}{\overset{\mathrm{\∞}}{\∫}}(Cx-r)dt$   式五

其中，r为系统输出参考量；显然有：

$\stackrel{\·}{\mathrm{\ρ}}=Cx-r$   式六

则可根据步骤2得到的最优解 $U_t^{*}={\[U_{t_0}{U}_{t_1}...U_{t_{N-1}}\]}^T,$ 求解其对应的状态量 ${\left[\begin{array}{c}x\\ \mathrm{\ρ}\end{array}\right]}_t^{*}=\left[\begin{array}{cccc}{x}_{t_0}& x_{t_1}& ...& x_{t_{N-1}}\\ {\mathrm{\ρ}}_{t_0}& {\mathrm{\ρ}}_{t_1}& ...& {\mathrm{\ρ}}_{t_{N-1}}\end{array}\right],$ 通过最小二乘法求取当前预测时域内T内最佳的比例、积分系数；

步骤4：在当前控制周期内，采用更新了控制系数的PI控制器实现对微网频率的调节控制。

2.根据权利要求1所述的一种基于MPC和PI控制方法的带电动汽车的微网调频方法，其特征在于，所述步骤2.2的具体方法是：根据预测控制的基本原理，首先以步骤2.1得到的最新测量值为初始条件，基于式三预测系统未来的动态；为此，设定预测时域为p，控制时域为m且m≤p；并定义：

定义一：控制时域之外，控制变量不变，即Δu(k+i)＝0,i＝m,m+1,…,p-1；

定义二：可测干扰在k时刻之后不变，即Δd(k+i)＝0,i＝1,2,…,p-1；

定义一是因为控制时域有可能小于预测时域，而预测系统未来动态需要在整个预测时域的控制输入；定义二是因为在当前k时刻，尚未获取干扰的未来取值；

在当前时刻k，测量值为x(k)，计算Δx(k)＝x(k)-x(k-1)；这个Δx(k)将作为预测系统未来动态的起点；由式三预测k+1到k+3时刻的状态增量如下：

Δx(k+1|k)＝AΔx(k)+B_uΔu(k)+B_dΔd(k)  式七

Δx(k+2|k)＝AΔx(k+1|k)+B_uΔu(k+1)+B_dΔd(k+1)  式八

＝A²Δx(k)+AB_uΔu(k)+B_uΔu(k+1)+AB_dΔd(k)

Δx(k+3|k)＝AΔx(k+2|k)+B_uΔu(k+2)+B_dΔd(k+2)  式九

＝A³Δx(k)+A²B_uΔu(k)+AB_uΔu(k+1)+B_uΔu(k+2)+A²B_dΔd(k)

上式中，k+1|k表示k时刻对k+1时刻的预测；符号|后面的k表示当前时刻为k；进而，可以预测k+m至k+p时刻的状态

Δx(k+m|k)＝AΔx(k+m-1|k)+B_uΔu(k+m-1)+B_dΔd(k+m-1)

＝A^mΔx(k)+A^m-1B_uΔu(k)+A^m-2B_uΔu(k+1)+…+B_uΔu(k+m-1)+A^m-1B_dΔd(k),

·

·  式十

·

Δx(k+p|k)＝AΔx(k+p-1|k)+B_uΔu(k+p-1)+B_dΔd(k+p-1)

＝A^pΔx(k)+A^p-1B_uΔu(k)+A^p-2B_uΔu(k+1)+…+A^p-mB_uΔu(k+m-1)+A^p-1B_dΔd(k),

进一步，有输出方程(2)预测k+1至k+p的被控输出

y_c(k+1|k)＝C_cΔx(k+1|k)+y_c(k)＝C_cAΔx(k)+C_cB_uΔu(k)+C_cB_dΔd(k)+y_c(k)  式十一

y_c(k+2|k)＝C_cΔx(k+2|k)+y_c(k+1|k)  式十二

＝(C_cA²+C_cA)Δx(k)+(C_cAB_u+C_cB_u)Δu(k)+C_cB_uΔu(k+1)+(C_cAB_d+C_cB_d)Δd(k)+y_c(k)

·

·

·

$\begin{array}{c}{y}_c(k+m|k)=C_c\mathrm{\Δ}x(k+m\left|k\right)+y_c(k+m-1|k)=\underset{i=1}{\overset{m}{\Σ}}{C}_c{A}^i\mathrm{\Δ}x\left(k\right)+\underset{i=1}{\overset{m}{\Σ}}{C}_c{A}^{i-1}\mathrm{\Δ}u\left(k\right)\\ +\underset{i=1}{\overset{m-1}{\Σ}}{C}_c{A}^{i-1}{B}_u\mathrm{\Δ}u(k+1)+\mathrm{...}+C_c{B}_u\mathrm{\Δ}u(k+m-1)+\underset{i=1}{\overset{m}{\Σ}}{C}_c{A}^{i-1}{B}_d\mathrm{\Δ}d\left(k\right)+y_c\left(k\right)\end{array},$   式十三

·

·

·

$\begin{array}{c}{y}_c(k+p|k)=C_c\mathrm{\Δ}x(k+p\left|k\right)+y_c(k+p-1|k)=\underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^i\mathrm{\Δ}x\left(k\right)+\underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^{i-1}{B}_u\mathrm{\Δ}u\left(k\right)+\underset{i=1}{\overset{p-1}{\Σ}}{C}_c{A}^{i-1}{B}_u\mathrm{\Δ}u(k+1)\\ +\mathrm{...}+\underset{i=1}{\overset{p-m+1}{\Σ}}{C}_c{A}^{i-1}{B}_u\mathrm{\Δ}u(k+m-1)+\underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^{i-1}{B}_d\mathrm{\Δ}d\left(k\right)+y_c\left(k\right)\end{array},$   式十四

定义p步预测输出向量和m输入向量如下：

$Y_p(k+1|k)\stackrel{def}{=}{\left[\begin{array}{c}{y}_c(k+1|k)\\ y_c(k+2|k)\\ .\\ .\\ .\\ y_c(k+p|k)\end{array}\right]}_{p\×1},\mathrm{\Δ}U\left(k\right)\stackrel{def}{=}{\left[\begin{array}{c}\mathrm{\Δ}u\left(k\right)\\ \mathrm{\Δ}u(k+1)\\ .\\ .\\ .\\ \mathrm{\Δ}u(k+m-1)\end{array}\right]}_{p\×1}$

那么，对系统未来p步预测的输出可以由下面的预测方程计算：

Y_p(k+1|k)＝S_xΔx(k)+Γy_c(k)+S_dΔd(k)+S_uΔU(k)  式十五

其中：

$S_x={\left[\begin{array}{c}{C}_{c}A\\ \underset{i=1}{\overset{2}{\Σ}}{C}_c{A}^i\\ .\\ .\\ .\\ \underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^i\end{array}\right]}_{p\×1},\mathrm{\Γ}={\left[\begin{array}{c}{I}_{n_c\×n_c}\\ I_{n_c\×n_c}\\ .\\ .\\ .\\ I_{n_c\×n_c}\end{array}\right]}_{p\×1},S_d{\left[\begin{array}{c}{C}_c{B}_d\\ \underset{i=1}{\overset{2}{\Σ}}{C}_c{A}^{i-1}{B}_d\\ .\\ .\\ .\\ \underset{i=1}{\overset{p}{\Σ}}{C}_c{A}^{i-1}{B}_d\end{array}\right]}_{p\×1},$