阿古斯蒂娜·斯基亚尔斯基 | 尼古拉斯·法埃多 | 约翰·V·林伍德
爱尔兰基尔代尔郡梅努斯大学电子工程系海洋能源研究中心
**摘要**
潮汐拦河坝发电厂利用海平面的日常变化来发电。优化潮汐拦河坝的运行需要解决一个相关的连续时间最优控制问题，这需要考虑过去的潮汐高度以及后续潮汐周期中预测的潮汐高度。尽管潮汐变化高度可预测，但它们也受到随机天气条件的影响，而在评估潮汐拦河坝的能源发电时，这些因素往往被忽视。本研究提出了一种递减视界控制框架，通过该框架可以计算潮汐拦河坝的最优运行方式，同时考虑过去的潮汐高度观测数据以及基于天气数据或仅基于天文潮汐的潮汐高度预测。在递减视界算法的每一步中，都使用基于矩的控制方法来计算潮汐拦河坝的最优运行，这种方法能够将无限维的潮汐拦河坝最优控制问题转化为有限维的非线性程序。通过这一框架，评估了在控制计算中考虑风暴潮的影响，与仅考虑天文潮汐相比，其效果有所提升，突显了在计算潮汐拦河坝最优运行时考虑风暴潮的价值。
**1. 引言**
部署可再生能源发电厂对于减少温室气体排放和可持续应对气候变化至关重要。在各种可再生能源中，潮汐范围（即高低海平面之间的差异）变化缓慢且高度可预测，因为大部分潮汐力来源于地球与月球以及地球与太阳之间的相互作用。潮汐范围的潜在能量可以通过潮汐拦河坝发电厂加以利用。图1展示了潮汐拦河坝的示意图，其中堤坝上装有涡轮机和水闸，将一个盆地与开阔海域分隔开来。随着潮汐水平的上升和下降，盆地内外的水位差驱动水流通过涡轮机，产生机械能，然后由与涡轮机相连的发电机将机械能转换为电能。沿拦河坝壁设置的水闸通过增加水流通过拦河坝的流量来提高运行的灵活性，而无需进行能量转换。此外，在某些时期，涡轮机还可以用于将水泵入或泵出盆地，从而增加可用的潜在能量和后续的发电量。大规模运行的潮汐拦河坝的例子包括法国的拉朗斯潮汐拦河坝和韩国的泗湖潮汐拦河坝；另一方面，英国沿海也提出了多个项目[1]，那里的潮汐范围是全球最高的之一。
在潮汐拦河坝中，涡轮机的发电量取决于涡轮机流量Qt和水头H（即盆地内水位ni与外部海平面no之间的差异）。ni取决于拦河坝之前的运行情况，即盆地之前被填充或排空的程度，从而影响H；此外，某一时刻的涡轮机流量Qt不仅决定了瞬时发电量，还决定了未来可用的水头H，即后续周期内可以产生的电量。因此，优化潮汐拦河坝的运行应被构建为一个时间窗口内的最优控制问题（OCP），其中需要包括过去和未来的潮汐周期，这表明需要潮汐高度预测。
**2. 潮汐拦河坝模型**
本节描述了用于控制计算的潮汐拦河坝模型。所研究的案例基于位于法国圣马洛的拉朗斯潮汐拦河坝。表1展示了本研究中考虑的拉朗斯潮汐拦河坝的设计参数[10]。
**3. 递减视界基于矩的控制框架**
本节介绍了递减视界基于矩的控制框架。第4节展示了使用不同潮汐预测获得的结果，第6节总结了结论。
**4. 结论**在基于矩的控制中，最优控制问题（OCP）可以用一组（三角函数）基函数来参数化，这些基函数由动态信号生成器生成，系统的稳态响应可以通过使用矩的数学概念来近似。如果问题被适当地提出，那么在矩域中近似的解保证会收敛到系统的稳态响应[15]。有兴趣的读者可以参考[2]、[16]以获取更多关于应用于潮汐坝的基于矩的控制的信息。
3.1. 潮汐坝的基于矩的控制
本文使用基于矩的控制框架将潮汐坝操作的无限维OCP转换为有限维的非线性程序，遵循在[2]中详细描述的过程。为了保持手稿的相对完整性，在本节中简要回顾了潮汐坝基于矩的参数化的数学表述。考虑描述单输入单输出（SISO）非线性潮汐坝系统的状态方程f:R×R×R²→R（通过将方程（4）代入方程（5））：(7)
ni̇=−Qt−sign(ni−no)CdsAs²g|ni−no|Ab(ni)=f(ni,no,u)。
利用主要由天文强迫驱动的潮汐高度的谐波行为，制定了以下谐波信号生成器G：(8)
G:θ̇=Sθ,no≈ne=Lneθ，
其中 (9)
θ=1cosωo×sinωo×...cosdωo×sinωo，
S=0⊕⨁p=1d₀pωo−pωo⁰，
θ∈Rⁿ，S∈Rⁿ×ν，ν=2d+1，ne是控制器内实际潮汐高度的描述。LneT∈Rⁿ是一个系数向量，它将ne投影到θ上；也就是说，潮汐高度no被近似为基频为ωo的谐波子空间的线性组合，d是谐波的数量。基频取决于控制计算所考虑的时间周期Th，即ωo=2π/Th；也就是说，问题的离散化取决于Th的选择，这将在第4.1节中讨论。对于一个固定的截止频率ωco，信号生成器G中存在的谐波数量为d=floor(ωco/ωo)。在潮汐坝最优控制的特定情况下，发现ωco应至少为20周期/天[2]；在本研究中，选择了ωco=28周期/天的值。
现在假设控制变量，由向量u=[Qt,As]表示，也可以用谐波函数来近似，假设天文潮汐决定了系统的主导动态。那么，控制变量被参数化为θ中基函数的线性组合：(10)
u=QtAs=LQtLAsθ=Luθ，
其中LuT=[LQtTLAsT]∈Rⁿײ是定义Qt和As在矩域中投影的向量。由此产生的互联系统具有以下形式 (11)
θ̇=Sθ,ni̇=f(ni,Lneθ,Luθ)。
假设信号生成器的初始条件θ(0)使得(Lne,Lu,S,θ(0))最小，并且考虑到S的特征值是纯虚数，存在一个映射π，即潮汐坝系统的矩，它解决了以下偏微分方程 (12)
∂π(θ)∂θSθ=f(π(θ),Lneθ,Luθ)。
此外，对于任何固定的轨迹θ(t)，系统的稳态响应可以根据其相关的矩来计算：(13)
niss(t)≈π(θ(t))=Lniθ(t),ni̇ss(t)≈Lniθ̇(t)=LniSθ(t)，
其中LniT∈Rⁿ是近似系统状态ni在矩域中的向量。因此，使用基于矩的参数化，连续时间的潮汐坝系统被离散化为由G给出的谐波函数的线性组合：(14)
θ̇=Sθ,no≈Lneθ,u=Luθ,ni≈Lniθ
Lni是通过应用类似Galerkin的方法[17]计算得出的，其中残差函数R在有限的时间配置点Nc∈N处被强制为0。R是通过将f(ni,ne,u)（方程（7）中的状态和输入）替换为它们的基于矩的参数化来定义的：(15)
R(Lni,Lne,Lu,t)=AbLniSθ(t)+LQtθ(t)+sign(Lni−Lne)θ(t)CdsLAsθ(t)²g|Lni−Lne|θ(t)。
强制R在Nc个配置点处为0，导致以下非线性方程组：(16)
R(Lni,Lne,Lu)=R(Lni,Lne,Lu,t₁)⋮R(Lni,Lne,Lu,tNc)=0，
只要Nc≥ν，就可以计算出Lni。
3.2. 递减视界控制
递减视界控制包括在时间间隔内迭代求解OCP，并为较短的控制周期检索控制变量的最优轨迹。本研究中实现的递减视界控制策略总结在以下伪代码算法中：
首先，在步骤N，外部输入ne在长度为Th的时间间隔ΦN内被检索，即时间视界，然后如第3.3节所述将其投影到信号生成器G上。ΦN以当前时间点tN为中心，即ΦN=ΦpN∪ΦfN，其中ΦpN=[tN−Th²,tN)是过去的时间间隔，ΦfN=[tN,tN+Th²]是未来的时间间隔。使用第3.1节中描述的基于矩的控制参数化，每个步骤N的潮汐坝OCP被如下制定：(17)
LuN,opt=argmax(LuN∫ΦNP(LniN,LneN,LuN)dt，
受限于
R(LniN,LneN,LuN)=0,
LniNA≤Bni,LuNA≤Bu,
ρgμN⊙(LniN−LnoN)Λ⊙LQtNΛ≤BPt，
其中μN是涡轮效率，在矩域参数化中映射（即作为LniN、LneN和LuN的函数），⊙表示标准的逐元素（Hadamard）乘积，(18)
Λ=θ(t₁)…θ(tNc),
A=ΛΛ,
Bni=nimax₁Nc/min₁Nc,
Bu=umax₁Nc/min₁Nc,
BPt=Ptmax₁Nc。
注意，在(17)中，约束在Nc个配置点处得到执行。然后为控制周期ΦcN=[tN,tN+Tc]（长度为Tc
两个连续控制解之间的不连续性通过在时间间隔φtn=(tn−1+tc,tn−1+tc+tt]=(tn,tn+tt]∈φcn内将步骤n−1的解与步骤n的解重叠来处理，长度为tt≪tc，并通过过渡函数[18]定义平滑过渡，如图4所示。结果控制解utn,opt计算如下：(19)
utn,opt=ξ(−t)un−1,optφtn+ξ(t)un,optφtn，
其中ξ:φtn→[0,1]是在t=t−tt/2处评估的双曲正切函数，定义为：(20)
ξ(t)=tanh(t)+1/2，
其中ξ(−t)+ξ(t)=1∀t。
3.3. 输入参数化
因为ocp是为以当前时间点tn为中心的时间周期φn求解的，输入ne可以使用实际观测值，即ne≔no，对于所有t∈φp，以及预测值，对于所有t∈φf，如图5所示。这样，关于t∈φp的实际潮汐变化的信息被包括在内，并在递减视界算法的每一步中更新。此外，递减视界框架允许使用t∈φf的潮汐高度预测，该预测也可以在每一步中更新。
为了确保基于矩的参数化中的ocp是适定的，系统的输入必须在用于控制计算的周期th内是周期性的。为此，在递减视界算法的每一步n中，ne被窗口化到一个th周期的映射ln：φn→[1,0]，该映射在极端处将ne平滑地驱动为零（见图5，右侧轴）。然后将窗口化的潮汐高度neln投影到由信号生成器g给出的谐波空间。因此，每个步骤n的neln的基于矩的参数化定义为：(21)
ne≈neln=lnne≈lnenθ。
图5显示了在φn的极端处，由于使用了窗口化映射，ne在矩域中的近似相对于实际潮汐高度no较差。另一方面，在φn的中心，ln对ne的表示的影响可以忽略不计。因此，用于控制计算的周期th影响了在φcn（对于固定的控制长度tc）问题的转录准确性，以及传递给控制器的关于外部输入的信息量（类似于标准的mpc方法）。此外，如第3.1节所述，信号生成器g的基频ωo取决于所选的th；th越长，对于固定的ωco，基于矩的参数化的维度越高，从而更准确地表示系统，但计算成本也更高。
图6显示了描述递减视界基于矩的控制器流程的框图。为了模拟外部输入ne，在预测视界φfn，考虑了三种情况，总结在表2中：
• 使用从st. malo的潮汐计观测中提取的观测潮汐no；
• 使用st. malo的eot20全球潮汐模型中的天文潮汐na（见第2.1节）；
• 使用在edf france网站上每周发布的st. malo的天气信息潮汐预测nf[19]。
4. 结果
本节展示了使用第3节中描述的递减视界基于矩的控制框架解决潮汐坝ocp所获得的结果。控制时间tc是从递减视界算法的每一步n检索控制解的时间，为1小时，nmax=720，因此模拟跨度为ω=1个月，如图7所示。然后如第3.2节所述，将控制解组合起来，创建一个最优轨迹uopt=[qtopt,asopt]，覆盖整个周期ω，然后在时间域中实现它，以计算系统的稳态响应ni和产生的能量e=∫ωptdt。
4.1. 控制视界长度th的选择
作为第一步，通过计算不同th值下的最优涡轮流量和水闸门面积来评估控制视界长度th对控制器性能的影响。在本节中，考虑了理想情况，即使用观测到的潮汐no来模拟过去（在φpn）和未来（在φfn）的潮汐高度。图8的左侧轴显示了使用每个th值计算出的最优运行产生的能量生成e；虚线标记e=97%emax，emax是在th=26小时时达到的最大能量，表明在th=24小时后，e保持在3%的范围内。请注意，本节显示的结果是针对选定的控制周期tc=1小时特定的。图8的右侧轴显示了状态ni与约束nimax和nimin的累积偏差δni，计算如下：(22)
δni=∫ωmax(ni−nimax,0)dt−∫ωmin(ni−nimin,0)dt。
由于基于矩的参数化是对系统的近似，因此预计在时间域中实现控制解（即qt和as的最优轨迹）时，系统状态ni的轨迹可能会偏离在控制器内计算的系统的矩lninθ。这种近似误差可能会通过递减视界算法的后续步骤传播，导致状态ni超出约束。较大的th值增强了在周期th中心的系统在矩域参数化中的近似，如第3.3节所述，减少了ni和lninθ之间的差异，从而减少了δni。图9显示了不同th值下系统状态的矩域参数化lninθ（黑色）和系统状态的时间域解ni（绿色），时间片段为420小时。随着th的增加，ni更接近lninθ，并且较少超出约束（由红色虚线表示）。预计e会随着th的增加而增加；如第3.3节所述，th的选择定义了传递给控制器的关于潮汐高度的（过去和未来的）信息量，以及ocp的基于矩的参数化的解空间。然而，图8显示e并不随th单调增加。这可以归因于控制解对每个th值导致的离散化的敏感性。la rance的潮汐高度中最主要的频率是对应于天文成分m2的频率，其周期为12.42小时，频率为ωm2=0.5059 rad/h。图10表示每个th值的解空间θ，其中每个点对应一个频率，每种颜色对应不同的th。对于每个解空间，特定频率与ωm2的距离最小（在图10中用红色边缘表示），计算为两个频率之间差异的l1范数：(23)
d(ω,ωm2)=|ω−ωm2|。
图10的右侧轴显示了d(ω,ωm2)随th的变化情况，类似于图8中e和δni的行为。距离d(ω,ωm2)在th=8小时时最大，此时e最低，δni最高。 两个连续控制解之间的不连续性通过在时间间隔φtn=(tN−1+Tc,tN−1+Tc+Tt]=(tN,tN+Tt]∈ΦcN内将步骤N−1的解与步骤N的解重叠来处理，长度为Tt≪Tc，并通过过渡函数[18]定义平滑过渡，如图4所示。结果控制解utN,opt计算如下：(19) utn,opt=ξ(−t)uN−1,optΦtN+ξ(t)uN,optΦtN， 其中ξ:φtn→[0,1]是在t=t−Tt/2处评估的双曲正切函数，定义为：(20) ξ(t)=tanh(t)+1/2， 其中ξ(−t)+ξ(t)=1∀t。 3.3. 输入参数化 因为ocp是为以当前时间点tn为中心的时间周期φn求解的，输入ne可以使用实际观测值，即ne≔no，对于所有t∈φp，以及预测值，对于所有t∈φf，如图5所示。这样，关于t∈φp的实际潮汐变化的信息被包括在内，并在递减视界算法的每一步中更新。此外，递减视界框架允许使用t∈φf的潮汐高度预测，该预测也可以在每一步中更新。 为了确保基于矩的参数化中的ocp是适定的，系统的输入必须在用于控制计算的周期th内是周期性的。为此，在递减视界算法的每一步n中，ne被窗口化到一个th周期的映射ln：φn→[1,0]，该映射在极端处将ne平滑地驱动为零（见图5，右侧轴）。然后将窗口化的潮汐高度neln投影到由信号生成器g给出的谐波空间。因此，每个步骤n的neln的基于矩的参数化定义为：(21) ne≈neln=LNne≈LneNθ。 图5显示了在φn的极端处，由于使用了窗口化映射，ne在矩域中的近似相对于实际潮汐高度no较差。另一方面，在φn的中心，ln对ne的表示的影响可以忽略不计。因此，用于控制计算的周期th影响了在φcn（对于固定的控制长度tc）问题的转录准确性，以及传递给控制器的关于外部输入的信息量（类似于标准的mpc方法）。此外，如第3.1节所述，信号生成器g的基频ωo取决于所选的th；th越长，对于固定的ωco，基于矩的参数化的维度越高，从而更准确地表示系统，但计算成本也更高。 图6显示了描述递减视界基于矩的控制器流程的框图。为了模拟外部输入ne，在预测视界φfn，考虑了三种情况，总结在表2中： • 使用从st. malo的潮汐计观测中提取的观测潮汐no； • 使用st. malo的eot20全球潮汐模型中的天文潮汐na（见第2.1节）； • 使用在edf france网站上每周发布的st. malo的天气信息潮汐预测nf[19]。 4. 结果 本节展示了使用第3节中描述的递减视界基于矩的控制框架解决潮汐坝ocp所获得的结果。控制时间tc是从递减视界算法的每一步n检索控制解的时间，为1小时，nmax=720，因此模拟跨度为Ω=1个月，如图7所示。然后如第3.2节所述，将控制解组合起来，创建一个最优轨迹uopt=[Qtopt,Asopt]，覆盖整个周期Ω，然后在时间域中实现它，以计算系统的稳态响应ni和产生的能量E=∫ΩPtdt。 4.1. 控制视界长度th的选择 作为第一步，通过计算不同th值下的最优涡轮流量和水闸门面积来评估控制视界长度th对控制器性能的影响。在本节中，考虑了理想情况，即使用观测到的潮汐no来模拟过去（在φpn）和未来（在φfn）的潮汐高度。图8的左侧轴显示了使用每个th值计算出的最优运行产生的能量生成e；虚线标记e=97%Emax，Emax是在Th=26小时时达到的最大能量，表明在Th=24小时后，E保持在3%的范围内。请注意，本节显示的结果是针对选定的控制周期Tc=1小时特定的。图8的右侧轴显示了状态ni与约束nimax和nimin的累积偏差Δni，计算如下：(22) δni=∫Ωmax(ni−nimax,0)dt−∫Ωmin(ni−nimin,0)dt。 由于基于矩的参数化是对系统的近似，因此预计在时间域中实现控制解（即qt和as的最优轨迹）时，系统状态ni的轨迹可能会偏离在控制器内计算的系统的矩lninθ。这种近似误差可能会通过递减视界算法的后续步骤传播，导致状态ni超出约束。较大的th值增强了在周期th中心的系统在矩域参数化中的近似，如第3.3节所述，减少了ni和lninθ之间的差异，从而减少了δni。图9显示了不同th值下系统状态的矩域参数化lninθ（黑色）和系统状态的时间域解ni（绿色），时间片段为420小时。随着th的增加，ni更接近lninθ，并且较少超出约束（由红色虚线表示）。预计e会随着th的增加而增加；如第3.3节所述，th的选择定义了传递给控制器的关于潮汐高度的（过去和未来的）信息量，以及ocp的基于矩的参数化的解空间。然而，图8显示e并不随th单调增加。这可以归因于控制解对每个th值导致的离散化的敏感性。la rance的潮汐高度中最主要的频率是对应于天文成分m2的频率，其周期为12.42小时，频率为ωm2=0.5059 rad h。图10表示每个th值的解空间θ，其中每个点对应一个频率，每种颜色对应不同的th。对于每个解空间，特定频率与ωm2的距离最小（在图10中用红色边缘表示），计算为两个频率之间差异的l1范数：(23) d(ω,ωm2)=|ω−ωM2|。 图10的右侧轴显示了d(ω,ωm2)随th的变化情况，类似于图8中e和δni的行为。距离d(ω,ωm2)在th=>两个连续控制解之间的不连续性通过在时间间隔φtn=(tn−1+tc,tn−1+tc+tt]=(tn,tn+tt]∈φcn内将步骤n−1的解与步骤n的解重叠来处理，长度为tt≪tc，并通过过渡函数[18]定义平滑过渡，如图4所示。结果控制解utn,opt计算如下：(19)
utn,opt=ξ(−t)un−1,optφtn+ξ(t)un,optφtn，
其中ξ:φtn→[0,1]是在t=t−tt/2处评估的双曲正切函数，定义为：(20)
ξ(t)=tanh(t)+1/2，
其中ξ(−t)+ξ(t)=1∀t。
3.3. 输入参数化
因为ocp是为以当前时间点tn为中心的时间周期φn求解的，输入ne可以使用实际观测值，即ne≔no，对于所有t∈φp，以及预测值，对于所有t∈φf，如图5所示。这样，关于t∈φp的实际潮汐变化的信息被包括在内，并在递减视界算法的每一步中更新。此外，递减视界框架允许使用t∈φf的潮汐高度预测，该预测也可以在每一步中更新。
为了确保基于矩的参数化中的ocp是适定的，系统的输入必须在用于控制计算的周期th内是周期性的。为此，在递减视界算法的每一步n中，ne被窗口化到一个th周期的映射ln：φn→[1,0]，该映射在极端处将ne平滑地驱动为零（见图5，右侧轴）。然后将窗口化的潮汐高度neln投影到由信号生成器g给出的谐波空间。因此，每个步骤n的neln的基于矩的参数化定义为：(21)
ne≈neln=lnne≈lnenθ。
图5显示了在φn的极端处，由于使用了窗口化映射，ne在矩域中的近似相对于实际潮汐高度no较差。另一方面，在φn的中心，ln对ne的表示的影响可以忽略不计。因此，用于控制计算的周期th影响了在φcn（对于固定的控制长度tc）问题的转录准确性，以及传递给控制器的关于外部输入的信息量（类似于标准的mpc方法）。此外，如第3.1节所述，信号生成器g的基频ωo取决于所选的th；th越长，对于固定的ωco，基于矩的参数化的维度越高，从而更准确地表示系统，但计算成本也更高。
图6显示了描述递减视界基于矩的控制器流程的框图。为了模拟外部输入ne，在预测视界φfn，考虑了三种情况，总结在表2中：
• 使用从st. malo的潮汐计观测中提取的观测潮汐no；
• 使用st. malo的eot20全球潮汐模型中的天文潮汐na（见第2.1节）；
• 使用在edf france网站上每周发布的st. malo的天气信息潮汐预测nf[19]。
4. 结果
本节展示了使用第3节中描述的递减视界基于矩的控制框架解决潮汐坝ocp所获得的结果。控制时间tc是从递减视界算法的每一步n检索控制解的时间，为1小时，nmax=720，因此模拟跨度为ω=1个月，如图7所示。然后如第3.2节所述，将控制解组合起来，创建一个最优轨迹uopt=[qtopt,asopt]，覆盖整个周期ω，然后在时间域中实现它，以计算系统的稳态响应ni和产生的能量e=∫ωptdt。
4.1. 控制视界长度th的选择
作为第一步，通过计算不同th值下的最优涡轮流量和水闸门面积来评估控制视界长度th对控制器性能的影响。在本节中，考虑了理想情况，即使用观测到的潮汐no来模拟过去（在φpn）和未来（在φfn）的潮汐高度。图8的左侧轴显示了使用每个th值计算出的最优运行产生的能量生成e；虚线标记e=97%emax，emax是在th=26小时时达到的最大能量，表明在th=24小时后，e保持在3%的范围内。请注意，本节显示的结果是针对选定的控制周期tc=1小时特定的。图8的右侧轴显示了状态ni与约束nimax和nimin的累积偏差δni，计算如下：(22)
δni=∫ωmax(ni−nimax,0)dt−∫ωmin(ni−nimin,0)dt。
由于基于矩的参数化是对系统的近似，因此预计在时间域中实现控制解（即qt和as的最优轨迹）时，系统状态ni的轨迹可能会偏离在控制器内计算的系统的矩lninθ。这种近似误差可能会通过递减视界算法的后续步骤传播，导致状态ni超出约束。较大的th值增强了在周期th中心的系统在矩域参数化中的近似，如第3.3节所述，减少了ni和lninθ之间的差异，从而减少了δni。图9显示了不同th值下系统状态的矩域参数化lninθ（黑色）和系统状态的时间域解ni（绿色），时间片段为420小时。随着th的增加，ni更接近lninθ，并且较少超出约束（由红色虚线表示）。预计e会随着th的增加而增加；如第3.3节所述，th的选择定义了传递给控制器的关于潮汐高度的（过去和未来的）信息量，以及ocp的基于矩的参数化的解空间。然而，图8显示e并不随th单调增加。这可以归因于控制解对每个th值导致的离散化的敏感性。la rance的潮汐高度中最主要的频率是对应于天文成分m2的频率，其周期为12.42小时，频率为ωm2=0.5059 rad/h。图10表示每个th值的解空间θ，其中每个点对应一个频率，每种颜色对应不同的th。对于每个解空间，特定频率与ωm2的距离最小（在图10中用红色边缘表示），计算为两个频率之间差异的l1范数：(23)
d(ω,ωm2)=|ω−ωm2|。
图10的右侧轴显示了d(ω,ωm2)随th的变化情况，类似于图8中e和δni的行为。距离d(ω,ωm2)在th=8小时时最大，此时e最低，δni最高。>随着Th的增加，d(ω,ωM2)减小，但在Th=14小时和Th=22小时之间有相对的增加；在相同的Th值下，E略有减小而Δni增加。在Th=26小时后，d(ω,ωM2)再次增加，但由于θ中的频率数量增加，E保持在Emax的3%范围内。下载：下载高分辨率图像（237KB）下载：下载全尺寸图像图8. 左轴：使用不同视界长度Th（蓝色方块）产生的能量E，以及与最大能量之间的3%范围（虚线）。右轴：盆地水平的累积约束偏差Δni。（关于此图例中颜色的解释，请参阅本文的网页版本。）下载：下载高分辨率图像（1004KB）下载：下载全尺寸图像图9. 系统状态的矩域参数化LniNθ（黑线），对于N∈[1,Nmax]，以及系统状态的时域解ni（绿线），对于不同的Th值。红色虚线代表约束nimax和nimin。（关于此图例中颜色的解释，请参阅本文的网页版本。）下载：下载高分辨率图像（210KB）下载：下载全尺寸图像图10. 左轴：基于矩的潮汐拦河坝系统的解空间θ，对于不同的Th值（用不同的颜色表示）。每个点代表θ中存在的一个频率；边缘为红色的点是最接近M2成分频率ωM2的频率。右轴：边缘为红色的频率与ωM2之间的距离。（关于此图例中颜色的解释，请参阅本文的网页版本。）4.2. 使用天文和天气信息预测在本节中，用于递归视界矩基控制器的输入ne结合了过去间隔ΦpN的观测值no和未来间隔ΦfN的预测值，对于所有N∈[1,Nmax]。考虑的两种预测是天文潮汐na和天气信息预测nf，如表2所示。图11用蓝色显示了在第4.1节中提到的在ΦfN使用no的理想情况，而使用na和nf作为预测的结果分别用绿色和黄色显示。图12显示了理想情况（Eno）与使用na作为预测（Ena）之间、理想情况与使用nf作为预测（Enf）之间以及Enf与Ena之间的能量百分比差异。理想情况与使用na之间的差异（图12(a)），解释为风暴潮对最优控制解决方案的影响，在1.5%到11.6%之间；另一方面，使用nf产生的能量更接近理想情况（图12(b)），与使用na相比；使用nf达到的最大能量比理想情况下的最大能量低1.4%，而使用na作为预测达到的最大能量不到理想情况下最大能量的97%。理想情况与使用nf之间的差异，解释为预测误差对最优控制解决方案的残余影响，低于5.4%，并且在Th=30小时时几乎可以忽略不计。请注意，无论使用na还是nf作为预测，与理想情况相比的能量差异都随着Th的增加而减小，能量生成也增加；换句话说，控制器的性能随着Th的增加而提高，这与使用的预测无关。图12(c)显示，使用nf作为预测可以将能量从1%增加到发电厂容量的6%。也就是说，在所考虑的研究案例中，使用天气信息预测有净收益。请注意，如第2.1节所述，nf与no之间的RMSE为0.15米，是最大潮汐范围13米的1.2%，而na的RMSE为0.35米，是最大潮汐范围的2.7%。也就是说，nf与no之间的RMSE比nf与no之间的RMSE低57%。另一方面，使用nf达到的最大能量与Emax之间的差距为1.4%，比使用na达到的最大能量与Emax之间的差距3.3%低57%。正如预期的那样，如果可用的天气信息预测的误差比天文潮汐大，那么控制器的性能可能会比简单地忽略风暴潮要差。因此，应使用na作为no的近似值的准确性作为评估天气信息预测是否适合控制计算的阈值。下载：下载高分辨率图像（250KB）下载：下载全尺寸图像图11. 考虑使用no、na和nf在ΦfN建模外部输入ne时的能量E（表2），使用不同的视界长度Th。（关于此图例中颜色的解释，请参阅本文的网页版本。）下载：下载高分辨率图像（89KB）下载：下载全尺寸图像(a)。在ΦfN使用no和na之间的能量百分比差异。下载：下载高分辨率图像（87KB）下载：下载全尺寸图像(b)。在ΦfN使用no和nf之间的能量百分比差异。下载：下载高分辨率图像（87KB）下载：下载全尺寸图像(c)。在ΦfN使用nf和na之间的能量百分比差异。图12. 比较了使用na（a）和nf（b）作为预测时产生的能量，以及使用na和nf作为预测时产生的能量（c）。Eno是在ΦfN使用no产生的能量，Ena是在ΦfN使用na产生的能量，Enf是在ΦfN使用nf产生的能量。图13显示了使用观测值作为预测（实线蓝线，如图9所示）以及使用天文潮汐na（绿线）和天气信息预测nf（黄线）作为预测时计算出的内盆地水位ni。使用天气信息预测nf的解决方案与仅使用天文潮汐na的解决方案相比，更接近理想情况。为了说明使用nf时控制器性能的提高，请注意，在图3(b)中，在退潮（低潮）时刻t=600小时，风暴潮相对较高；因此，如图13所示，使用na作为预测时，对于所有Th值，内盆地水位ni都超过了最大约束，而使用nf得到的控制解决方案使ni保持在约束范围内。关于u=[Qt,As]的最优轨迹，可以看出，使用nf作为预测的解决方案[Qtnf,Asnf]在最小二乘意义上更接近使用观测值作为预测的解决方案[Qtno,Asno]，而不是使用na作为预测的解决方案[Qtna,Asna]。图14显示了理想情况下的最优解与每种预测之间的RMSE；线条代表涡轮流量Qt的RMSE（左轴），而方块代表泄洪闸门面积As的RMSE（右轴），对于不同的Th值。正如预期的那样，对于每个Th值，使用nf作为预测得到的最优控制轨迹的RMSE都低于使用na作为预测得到的轨迹。下载：下载高分辨率图像（793KB）下载：下载全尺寸图像图13. 使用观测到的潮汐no（理想情况）、天文潮汐na和天气信息预测nf作为预测时系统状态ni的时域解，对于不同的Th值。红色虚线代表约束nimax和nimin。（关于此图例中颜色的解释，请参阅本文的网页版本。）下载：下载高分辨率图像（234KB）下载：下载全尺寸图像图14. 左轴：使用na作为预测的控制解决方案Qtna（蓝色线）与理想控制解决方案Qtno（蓝色线）之间的RMSE，以及使用nf作为预测的控制解决方案Qtnf（黄色虚线）与Qtno（蓝色线）之间的RMSE。右轴：使用na作为预测的控制解决方案Asna（蓝色方块）与理想控制解决方案Asno（蓝色方块）之间的RMSE，以及使用nf作为预测的控制解决方案Asnf（黄色方块）与Asno（黄色方块）之间的RMSE。（关于此图例中颜色的解释，请参阅本文的网页版本。）4.3. 计算考虑和实时实现如第3节所述，增加视界Th会降低基本频率ωo；因此，对于固定的截止频率ωco，与矩域参数化相关的信号生成器中的维度会增加。从数值上讲，这种频率的增加意味着优化变量的数量增加，即控制向量Lu的系数数量2ν（每个控制向量LQt和LAs都有ν个系数）。此外，更高的维度需要更多的配置点Nc来求解残差函数（16）。为了满足递归视界控制器的实时实现需求，增加的计算运行时间不应超过Tc，这里假设为1小时。在本研究中，使用Matlab®进行了模拟，采用活动集求解器[20]，时间步长dt=1.2分钟。使用的PC配备了第13代Intel®CoreTM i7-1365U处理器和16 GB的RAM。表3显示了对于不同的Th值，优化变量2ν的数量、时间步长Nt（即Th/dt）以及配置点Nc的数量，Nc大约选择为2ν。图15显示了本节考虑的不同Th值下，递归视界控制的720次迭代的平均运行时间和最大运行时间（使用no、na和nf作为预测）。模拟运行时间使用Matlab®函数tic和toc进行测量。所有三种情况的平均模拟运行时间几乎随Th线性增加，并且即使Th=30小时也相对较低。最大运行时间始终低于50分钟，这对于实时实现是可以接受的；然而，如果天气信息预测在每个递归视界步骤都更新，这样的余量可能不足以满足所使用的天气模型的模拟时间。下载：下载高分辨率图像（344KB）下载：下载全尺寸图像图15. 不同Th下的平均运行时间（左图）和最大运行时间（右图），使用no、na和nf。表3. 不同Th的模拟参数。Th [h]810121416182022242628302ν38465866748694102114122130142Nt401501601701801901100111011201130114011501Nc394658677386961011151241281435. 扩展结果：英国卡迪夫泻湖项目为了展示所提出的控制框架在各种设置下的一致性，本节提供了另一个案例研究，即在英国塞文河口，那里的潮汐范围资源是世界上最丰富的之一[21]。由于潮汐范围资源和风暴潮都具有很强的地点特异性，因此这里的目的是在不同的位置测试控制器，并评估从La Rance研究案例获得的先前结果的普遍性。本节考虑的潮汐拦河坝发电厂基于在[22]中详细描述的卡迪夫泻湖项目；用于建模卡迪夫泻湖的设计参数总结在表4中。图16显示了卡迪夫的风暴潮分布，在-1.5米到1.5米之间。图17显示了在所研究的时间内，所提出的卡迪夫泻湖位置观测到的潮汐no、天文潮汐na和风暴潮ns。对于这里考虑的720小时周期，no和na之间的RMSE为0.48米，最大潮汐范围约为12米，最大绝对涌浪为1.3米，即大约是潮汐范围的11%（大于第2.1节中La Rance数据中的8%）。表4. 卡迪夫泻湖项目的参数化，来自[22]。参数值单位最大泄洪闸门面积As7200m²涡轮机数量nt98–涡轮机额定功率20MW最大涡轮机流量400m³/s泄洪闸门流量系数Cds1–潮汐拦河坝OCP使用与La Rance案例中相同的递归视界矩基控制框架进行求解，Tc=1小时，以及相同的Th值范围。因此，在时域参数化中使用的解空间对于La Rance案例和Cardiff Lagoon案例在每个Th值上都是相同的。如图18所示，使用实际观测数据（而不仅仅是天文潮汐na）作为预测的理想情况下，产生的能量要显著多于仅使用天文潮汐na作为预测的情况。对于每个考虑的Th值，使用na作为预测的结果产生的能量比理想情况下低5%到78%，所有能量值都低于理想情况下最大能量的93%。显然，不仅风暴潮的存在会影响拦河坝的最优控制解决方案（如La Rance案例所示），而且相对于潮汐范围而言，更高的风暴潮水平会导致仅使用天文潮汐作为预测时能量损失更大。
下载：下载高分辨率图像（137KB）
下载：下载全尺寸图像
图16. 圣马洛地区风暴潮驱动的潮位分布。数据来自2020年至2024年。
下载：下载高分辨率图像（280KB）
下载：下载全尺寸图像
图17. 卡迪夫泻湖位置观测到的潮位、作为控制器外部输入的天文潮汐预测以及风暴潮。
[9]中的研究调查了风暴潮对英国西海岸不同拟议项目潮汐能发电的影响，包括塞文河口。[9]的结果表明，在固定（未优化）的操作策略下，使用天文潮位作为输入的模拟能量比包括风暴潮时的能量低3%，尤其是在塞文河口的泻湖中。也就是说，风暴潮对能量发电的影响远低于本节所见的情况，其中提出的卡迪夫泻湖在使用天文潮位作为预测时最多只能产生理想情况下最大能量的93%。然而，[9]的结果与本研究不可比较，因为[9]中没有进行优化。
此外，图18显示，在18小时到24小时之间，能量有相对下降的趋势，这与La Rance案例（图8）的情况类似。在塞文河口，与圣马洛一样，M2成分是最主要的成分；因此，如第4.1节所解释的，当Th在16小时到22小时之间时，解空间θ无法很好地表示系统动态，因为d(ω,ωM2)（方程（23））增加。
下载：下载高分辨率图像（243KB）
下载：下载全尺寸图像
图18. 左轴：在卡迪夫泻湖案例中，不使用na和使用na在ΦfN时模拟的能量E，不同的预测范围Th。右轴：在ΦfN时使用na和不使用na之间的能量差异ΔE百分比。
6. 结论
本文提出了一种新颖的非线性滚动时域基于矩的控制框架，可以在处理潮汐拦河坝最优控制（OCP）问题时考虑天气驱动的潮汐变化。可以看出，预测范围Th显著影响控制解决方案，因为它直接关系到时域中输入和解空间的结构。鉴于M2成分在潮位中的主导地位，当问题的离散化包括接近M2成分频率的频率时，控制器的性能会得到改善。
仅使用天文潮汐作为预测，在滚动时域控制器中产生的能量低于具有完美潮位预测的理想情况。这表明，在所考虑的位置，尽管潮汐通常被认为是确定性的，但风暴潮对最优控制解决方案和随后的能量产生有不可忽视的影响。此外，结果表明，在研究的时间段内，当在潮位预测中加入天气驱动的潮汐变化预测时，与仅使用天文潮位相比，能量产生有所增加。这些结果强调了如果可能的话，将天气信息丰富的预测作为潮汐拦河坝OCP的输入的价值，前提是该预测在最小二乘意义上比天文潮汐更接近实际潮位。从这个意义上说，所提出的控制框架适用于潮汐拦河坝的在线调度，因为它可以有效地包含实时潮位信息以及更新的天气驱动风暴潮预测。滚动时域方法使得潮汐拦河坝的调度更加灵活，可以根据变化且往往不可预测的天气条件调整操作。此外，可以使用基于矩的控制框架高效计算最优拦河坝操作，从而实现每小时更新控制解决方案。
本研究的局限性之一是使用的天气驱动潮汐预测在滚动时域控制器的每次迭代中都不会更新。一般来说，天气预测的拟合度会随时间下降，这可能导致随着Th的增加而降低控制器性能。需要进一步的研究来分析使用定期更新预测时的结果，而这在本研究中并未实现。
cRediT作者贡献声明：
Agustina Skiarski：撰写——原始草稿、可视化、方法论、调查、形式分析、数据整理、概念化。
Nicolás Faedo：撰写——审阅与编辑、监督、形式分析。
John V. Ringwood：撰写——审阅与编辑、监督、概念化。
