# 美赛

## D题

可以参考的资料：

- [Plan 2014](https://ijc.org/sites/default/files/IJC_LOSR_EN_Web.pdf)
- 美赛 2022B
- [GLERL Great Lakes Monthly Hydrologic Data (1860-Recent)](https://www.glerl.noaa.gov/ahps/mnth-hydro.html)
- 这里有各个湖的降水、蒸发、流入量。
- https://www.glerl.noaa.gov/ftp/publications/tech_reports/glerl-109/

### 第一问






确定五大湖在一年中任意时间，考虑各种利益相关者的利益。每个利益相关者的重要程度不同。

首先需要找到哪些是利益相关方，随后需要对其进行量化。
每个利益相关者，对于水位的评价函数 P(H) 如何制定？可以考虑模糊综合评价里的那个指派法。


不同利益相关者之间的重要程度如何确定？若利益相关者较多，可以采用主成分分析法(PCA) 进行降维，使用往年的数据得到相应的满意度样本矩阵，随后提取出主成份。再进入下面的步骤处理。

若利益相关者较少，可以直接使用层次分析法来评判重要性。

随后即可得到 (评价函数，重要程度) 的二元组，作为我们的利益相关方的参考。

接下来就是决定最佳水位，跑一个单目标规划，自变量设置为水位，最大化满意度即可。

主要的难点在于：确定利益相关者的评价函数。

主要的六大 interests:

- Domestic, Municipal and Industrial Water Uses
- Commercial Navigation
- Hydroelectric Generation
- Ecosystems
- Costal Zone
- Recreational Boating and Tourism

忽略后俩，考察前四个主要的利益。

#### Domestic, Municipal and Industrial Water Uses

> Changing lake levels may impact each water withdrawal
facility differently, depending, among other factors, on the
location of the facility, the infrastructure of the intake, and
the amount of water withdrawn.
> High water conditions in the past have not been a
significant problem for water withdrawal facilities in the
upper Great Lakes basin. Some facilities have experienced
problems with flooding of buildings, tunnels and property,
as well as shoreline damages. High water conditions also
can affect wastewater treatment plants by increasing the
infiltration into the plant, thus increasing the demand for
water and temporarily increasing treatment costs.
> Lower lake levels can lead to insufficient water depths at
intakes, water quality and navigation problems at water
withdrawal facilities, damage to equipment, more problems
with algae at intakes, and wells not being deep enough
during high demand periods. These problems would be
made more difficult by rising water temperatures in the
changing climate.

看起来大概是，高了会导致污水处理之类的成本变高，低了会导致抽不上水。水位过低是更重大的危险，水位高一点是还可以接受的损失。

取往年均值记为 $avg$，我们使用 $\Gamma$ 型函数进行指派：

$$\mu _A = \begin{cases}
e^{k_1(x-a)}, x < a\\
1, a \le x \le b\\
e^{-k_2(x-b)}, x > b
\end{cases}$$

取 $a = 0.999avg$，$b = 1.002avg$，容忍一定程度的高水位，


假设 $avg=174$，在 $174.5$ 左右水位时，高水位导致利益减半。相应的，$173.7$ 左右的低水位就会导致利益减半。当出现较大波动时，超出 $1$ 米的水位差，在两侧都会导致利益的极速衰减。

左侧 $k=6$，右侧 $g=-4$，对应不同的敏感性。

对五个湖都套用这个函数，得到的评分做分析结果如下：

                 0
count  1260.000000
mean      0.915066
std       0.201266
min       0.046872
25%       1.000000
50%       1.000000
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.642141
std       0.368385
min       0.014742
25%       0.262615
50%       0.732545
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.680285
std       0.360878
min       0.002700
25%       0.337621
50%       0.909682
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.709407
std       0.347295
min       0.007311
25%       0.386371
50%       0.969516
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.563867
std       0.341738
min       0.003281
25%       0.261980
50%       0.543464
75%       0.980702
max       1.000000

#### Commercial Navigation

整体上会愿意大，但不能太大。

$$\mu _A = \begin{cases}
e^{k_1(x-a)}, x < a\\
1, a \le x \le b\\
e^{-k_2(x-b)}, x > b
\end{cases}$$

依然使用 $\Gamma$ 型函数。选择 $a = 0.999avg$，$b = 1.01avg$，$k_1 = 9$，$k_2 = -2$，对水位上升不敏感、下降敏感。


试评估结果如下：

                 0
count  1260.000000
mean      0.899789
std       0.240719
min       0.010148
25%       1.000000
50%       1.000000
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.710631
std       0.407932
min       0.001790
25%       0.225970
50%       1.000000
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.742027
std       0.390722
min       0.000140
25%       0.352752
50%       1.000000
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.759528
std       0.376814
min       0.000625
25%       0.407527
50%       1.000000
75%       1.000000
max       1.000000
                 0
count  1260.000000
mean      0.682536
std       0.396201
min       0.000188
25%       0.251772
50%       1.000000
75%       1.000000
max       1.000000

#### Hydroelectric Generation Interests

我们认为在一定限度内，水越多、水电站产生的利益越大。但是超过某个限度后开始下降。这个 turning point 取 $1.01 \times avg$.

$a=b=1.01 \times avg$，同时我们认为水电站相对而言不那么敏感。$k_1 = 0.1$，$k_2 = -0.2$.


评估结果：

                 0
count  1260.000000
mean      0.832596
std       0.016891
min       0.776654
25%       0.821388
50%       0.832968
75%       0.844711
max       0.874800
                 0
count  1260.000000
mean      0.838943
std       0.034373
min       0.767683
25%       0.810075
50%       0.838720
75%       0.864695
max       0.931109
                 0
count  1260.000000
mean      0.840083
std       0.033288
min       0.747424
25%       0.815362
50%       0.840194
75%       0.864052
max       0.927630
                 0
count  1260.000000
mean      0.840724
std       0.031042
min       0.760663
25%       0.817450
50%       0.840662
75%       0.862806
max       0.925367
                 0
count  1260.000000
mean      0.927634
std       0.030296
min       0.837303
25%       0.907039
50%       0.926288
75%       0.949738
max       0.999140

#### Ecosystem Interests

人为的管控会影响生态系统多样性，因为很多生物本身依赖周期性的水位变化。

也就是不喜欢平均(

定义历史最大偏差 $D = \max \operatorname{abs}(X - avg)$，某数据偏差值 $x$ 与历史最大偏差的比值 $\dfrac{x}{D}$ 可以用于衡量偏差程度，这是一个可以超过 $1$ 的量。我们认为一定程度的偏差是有益的，但依然会超出后有害。

$$1 - 0.5 |\ln (\dfrac{x}{D} + 0.2)|$$


                 0
count  1260.000000
mean      0.543185
std       0.195416
min       0.207631
25%       0.383682
50%       0.550388
75%       0.688812
max       0.993771
                 0
count  1260.000000
mean      0.637328
std       0.207316
min       0.196959
25%       0.479287
50%       0.666652
75%       0.797536
max       0.999933
                 0
count  1260.000000
mean      0.590387
std       0.203769
min       0.197038
25%       0.432324
50%       0.601333
75%       0.742061
max       0.999772
                 0
count  1260.000000
mean      0.611336
std       0.207149
min       0.205027
25%       0.452789
50%       0.620054
75%       0.774599
max       0.999519
                 0
count  1260.000000
mean      0.548183
std       0.194548
min       0.199606
25%       0.394771
50%       0.541437
75%       0.684334
max       0.997555

至此，我们的评价函数已经完成了。

#### 权值

熵权法结果不是很好，使用 AHP 层次分析法确定权值。



给出矩阵，A B C D 代表 domestic 家庭用水、navigation 航行、hydroelectric 水电、ecology 生态

1 3 4 5
1/3 1 4 5
1/4 1/4 1 1
1/5 1/5 1 1

一致性检验指标 $CI = 0.053579956780095195` 通过了一致性检验。

计算得出的权值：

[0.51289291]
[0.30707541]
[0.09492501]
[0.08510667]

单目标优化模型，湖的水位作为自变量，目标函数即为加权的总评价函数。认为每个湖享有平等的权利。

使用 scipy minimize 中的 `nelder-merd` 算法直接计算，得到最优水位为：

[183.78029519 176.66465083 175.39084977 174.52226052 74.9175519 ]

对比平均水位：

[183.41346825396826, 176.44929365079363, 175.04081746031744, 174.1739126984127, 74.76801587301588]

略有差异，看出似乎全都是往上了一点点(


问题不大。

### 第二问

到了这里就需要建立一个湖水的模型了。


虽然说是五大湖，但是我们得到的数据里面，有下面这些地方的水位：

- St. Mary's River - Flow
- Lake Michigan and Lake Huron - Mean Water Level
- St. Clair River - Flow
- Lake St. Clair - Mean Water Level
- Detroit River - Flow
- Lake Erie - Mean Water Level
- Niagara River - Flow at Buffalo
- Lake Ontario - Mean Water Level
- Ottawa River - Flow at Carillon
- St. Lawrence River - Flow at Cornwall

观察一下，可以发现，五大湖给出的都是 Mean Water Level 平均水位，而河流给出的都是 Flow.

Flow 的单位是 $\mathrm{m}^3 / \mathrm{s}$，而水位单位是 $\mathrm{m}$. 中间还有一个体积到水位高度的换算，每个湖都不太一样。

我们的两个关键控制点位：the Compensating Works，在 St.Marys River 那里。
Marie and the Moses-Saunders Dam，在 ST.LAWRENCE RIVER 那里。

从左到右，湖、河分别是：

- LAKE SUPERIOR
- **St. Marys River(Controlled by the Compensating Works)**
- LAKE HURON
- LAKE MICHIGAN(这两个湖比较特殊，需要看俯视图才能理解)
- LAKE HURON 流出至：
- St. Clair River
- Lake St.Clair
- Detroit River
- LAKE ERIE
- Niagara River
- LAKE ONTARIO
- **ST. LAWRENCE RIVER(Controlled by Marie and the Moses-Saunders Dam)**
- Atlantic Ocean

看给出的数据，Lake Michigan and Lake Huron 直接合并在一起给了，那我们直接合并一起处理就行了。


对于每个湖建模，湖下一时刻的水位状态由上一时刻的如下影响因素决定：

- 水位(蓄水量)
- 流入量
- 流出量

事实上，水位只是表象，蓄水量才是根本。那么我们还是先计算蓄水量。

每个湖维护一个蓄水量到水位的转换映射。

蓄水量由两个因素直接决定：

- 流入量
- 流出量

其中，流入量由以下因素决定：

- 降雨量
- 河流流入量
- Winter Snowpack 冰雪融雪补给
- Ice Jams 冰塞

流出量由以下因素决定：

- 蒸发量
- 河流流出量
- 人类用水量

其他因素以后再考虑吧(

我们可以控制的是两个大坝，大坝具体能做出什么样的调控呢？

控制水的流速(

大坝有如下的属性：

- 最大蓄水量(可以暂不考虑)
- 最大流速(全开的流速)

理论上关的时候可以全关，全开可以最大流速。

- [Soo Compensating Works Dam](http://superior.uslakes.info/DamInfo.asp?DamID=B50CBBD2-D7C6-473A-B7E2-1C0C8DE5D5BB)
- 最大流速：$1925 \mathrm{m}^3/\mathrm{s}$
- 最大蓄水量：$1.2 \times 10^{13} \mathrm{m}^3$
- Marie and the Moses-Saunders Dam
- 看到一个[新闻](https://ijc.org/en/loslrb/lake-ontario-outflows-be-maintained-above-navigation-limit)，最大流速可以取 $12000 \mathrm{m}^3/\mathrm{s}$

对于两个湖，如果不加任何调控，影响他们之间的河流流速的变量应该是：

- 时间(冰塞)
- 两个湖的高度差

先测测数据再说吧。

#### 河流流量(Discharge) 估算

可调控的河流我们自己决定流量。

不可调控的需要根据两侧的湖进行计算。

##### St. Clair River

上接 LAKE HURON & LAKE MICHIGAN，下接 Lake St.Clair.

根据这个链接查到了每天的河流量(嗯对，每天、、、)：https://waterdata.usgs.gov/nwis/dv?cb_30208=on&format=html&site_no=04159130&legacy=&referred_module=sw&period=&begin_date=2008-01-01&end_date=2024-02-01

然后是 Lake Huron & Michigan 的水位数据。

https://www.lre.usace.army.mil/Missions/Great-Lakes-Information/Great-Lakes-Information-2/Water-Level-Data/

在这里下载到了所有的水位数据，按月。所以上面的流量也按月了。

数据处理，直接月份、高度为自变量，流量为因变量，得到的相关系数矩阵如图：


感觉不是很好，分季节之后：

春季的相关系数矩阵：

              month    height  discharge
month      1.000000 -0.065565   0.148549
height    -0.065565  1.000000   0.745297
discharge  0.148549  0.745297   1.000000

春季的线性回归结果：


夏季的相关系数矩阵：

              month    height  discharge
month      1.000000  0.147961   0.030745
height     0.147961  1.000000   0.925719
discharge  0.030745  0.925719   1.000000

夏季的线性回归结果：


秋季的相关系数矩阵：

              month    height  discharge
month      1.000000  0.219979  -0.053460
height     0.219979  1.000000   0.885361
discharge -0.053460  0.885361   1.000000

秋季的线性回归结果：


冬季的相关系数矩阵：

              month    height  discharge
month      1.000000 -0.079939   0.324969
height    -0.079939  1.000000   0.318698
discharge  0.324969  0.318698   1.000000

冬季的线性回归结果：


总结：夏秋的拟合优度不错，春季还可以接受，冬季的拟合优度差。

不管怎样，先跳过这个奇怪的冬季。

##### Detroit River

这里可以找到河流量的数据：https://waterdata.usgs.gov/nwis/uv?site_no=04165710&legacy=1

湖的数据就用赛题里提供的了。

这个的拟合情况比上一条河还不如。。。标准差在 10% 往上了。

##### Niagara River

找数据，又是这个网站。[又是这个网站。](https://waterdata.usgs.gov/nwis/dv/?site_no=04216000&PARAmeter_cd=00060)

更加糟糕的拟合情况，只能说存在一定的相关性。。。

#### 主要建模

我们的调控自变量只有两个：St. Marys River 的流量和 ST. LAWRENCE RIVER 的流量。

我们的湖水量由以下因素决定：

- 基础量(上一时刻的湖水量) $base$
- 流入量
- 上一条河的流量 $flow_{in}$
- 降雨量 $perception$
- 流出量
- 下一条河的流量 $flow_{out}$
- 蒸发量 $evaporation$
- 其他环境因素(人类取水量等) $other$

### Midlakes 模型

NOAA 的公开论文 A Coordinated Hydrologic Response Model for the Middle Great Lakes.

1998 年提出的论文，带有 fortran 代码但是跑不起来，自己重写一份并做一些优化调教。

Quinn 1978 年提出的论文中，给出了三种解法，Midlakes 模型选择其中的 'hydrologic' finite difference solution.

第一个等式，outflow equation:

$$QO_{j,m} = QO_{j,t} + \dfrac{1}{2} \left( \dfrac{\partial QO_{j}}{\partial t}\right) \mathrm{d} t$$

其中 $QO$ 代表 outflow, $m$ 代表该 timestep 中的均值，$t$ 代表在时间 $t$ 时的 outflow, $j=0,1,2,3,4$ 代表 Lakes Superior, Michigan-Huron, St. Clair, Erie, and Upper Niagara River levels.

接下来是 stage-fall-discharge equation：

$$QO_{j,t} = K_j \left[ \varphi Z_{j,t} + (1-\varphi_j)Z_{j+1,t} - ym_j \right]^{a_j} \left( Z_{j,t} - Z_{j+1,t} \right)^{b_j}$$

其中 $Z$ 为 lake surface elevation, 也就是 water level.

$K$ 是一个经验性的 parameter，outflow equation coefficient.

$ym$ 是 mean channel bottom elevation, $a$,$b$ 是 depth、fall exponents.

$\varphi$ 是一个加权参数，允许在决定 outflow 时给上游、下游更大的权重。

可以发现，$K$、$ym$、$a$、$b$、$\varphi$ 都是未知的参数，需要跑最优化进行拟合。

最后一个式子：

$$P_{j,m} + R_{j,m} + QO_{j-1,m} - QR_{j-1,m} \pm D_m \pm G_{j,m} = EV_{j,m} + CU_{j,m} + QO_{j,m} - QR_{j,m} + A_j \left( \dfrac{\Delta Z ^{t \to t + \Delta t}_{j}}{\Delta t} \right)$$

P for precipitation, R for runoff, QR for ice-weed retardation(冰塞),EV for evaporation rate, D for diversion, G for groundwater, CU for consumptive use.

代入三个式子消元。


Lake Erie 有:

$$\begin{aligned}
& P_3 + R_3 - EV_3 \pm D_3 \pm G_3 - CU_3 + K_2[\varphi_2 Z_2 + (1-\varphi_2)Z_3 -ym_2]^{a_2}(Z_2 - Z_3)^{b_2} - \\
& QR_2 - K_3[\varphi_3Z_3 +(1-\varphi_3)Z_4 - ym_3]^{a_3}(Z_3-Z_4)^{b_3} + \\
& QR_3 + \{ \dfrac{1}{2} K_3b_3[\varphi_3Z_3 +(1-\varphi_3)Z_4 - ym_3]^{a_3}(Z_3-Z_4)^{b_3 - 1} - \\
& \dfrac{1}{2} K_3a_3(1-\varphi_3)[\varphi_3Z_3 +(1-\varphi_3)Z_4 - ym_3]^{a_3 - 1}(Z_3-Z_4)^{b_3} \} \Delta Z_4 \\
& = \\
& (0)\Delta Z_1 - \\
& \{\dfrac{1}{2} K_2b_2[\varphi_2 Z_2 + (1-\varphi_2)]Z_3 - ym_2]^{a_2}(Z_2 - Z_3)^{b_2 - 1} + \dfrac{1}{2} K_2a_2\varphi_2[\varphi_2 Z_2 + (1-\varphi_2)]Z_3 - ym_2]^{a_2-1}(Z_2 - Z_3)^{b_2} \} \Delta Z_2 + \\
& \{ \dfrac{A_3}{\Delta t} + \dfrac{1}{2} K_3b_3[\varphi_3 Z_3 + (1-\varphi_3)]Z_4 - ym_3]^{a_3}(Z_3 - Z_4)^{b_3 - 1} + \dfrac{1}{2} K_3a_3\varphi_3[\varphi_3 Z_3 + (1-\varphi_3)]Z_3 - ym_3]^{a_3-1}(Z_3 - Z_4)^{b_3} + \\
& \dfrac{1}{2} K_2b_2[\varphi_2 Z_2 + (1-\varphi_2)]Z_3 - ym_2]^{a_2}(Z_2 - Z_3)^{b_2 - 1} - \dfrac{1}{2} K_2a_2\varphi_2[\varphi_2 Z_2 + (1-\varphi_2)]Z_3 - ym_2]^{a_2-1}(Z_2 - Z_3)^{b_2}
\} \Delta Z_3
\end{aligned}$$


Lake Stclair 有


Lake Huron&Michigan 有


这一团团的式子有点过度复杂了，让我们尝试做一些简化。

$$N(i) = P_i + R_i \pm D_i \pm G_i - EV_i - CU_i$$

$$OF(i) = K_i[\varphi_i Z_i + (1-\varphi_i)Z_{i+1} - ym_i]^{a_i}(Z_i - Z_{i+1})^{b_i}$$

$$OFB(i) = \dfrac{1}{2}K_ib_i[\varphi_i Z_i + (1-\varphi_i)Z_{i+1} - ym_i]^{a_i}(Z_i - Z_{i+1})^{b_i-1}$$

$$OFA(i, t) = \dfrac{1}{2}K_it[\varphi_i Z_i + (1-\varphi_i)Z_{i+1} - ym_i]^{a_i-1}(Z_i - Z_{i+1})^{b_i}$$

常见的用法包括 $OFA(i, \varphi)$ 与 $OFA(i,1-\varphi)$.

重写我们的 Lake Erie:

$$\begin{aligned}N(3) + OF(2) - QR_2 - OF(3) + QR_3 + [OFB(3) - OFA(3, 1- \varphi_3)]\Delta Z_4 = \\
(0)\Delta Z_1 - [OFB(2) + OFA(2,\varphi_2)]\Delta Z_2 + \\
[\dfrac{A_3}{\Delta t} + OFB(3) + OFA(3,\varphi_3) + OFB(2) - OFA(2, 1-\varphi_2)]\Delta Z_3
\end{aligned}$$

Lake stclair 有：

$$\begin{aligned}
N(2) + OF(1) - QR_1 - OF(2) + QR_2 = \\
-[OFB(1) + OFA(1,\varphi_1)] \Delta Z_1 \\
+[\dfrac{A_2}{\Delta t} + OFB(2) + OFA(2,\varphi_2) - OFB(1) - OFA(1,1-\varphi_1)]\Delta Z_2\\
-[OFB(2)-OFA(2,1-\varphi_2)]\Delta Z_3
\end{aligned}$$

Lake Huron & Michigan 有：

$$\begin{aligned}
N(1) + QO_{0,m} - OF(1) + QR_1 = \\
[\dfrac{A_1}{\Delta t} + OFA(1,\varphi_1) + OFB(1)] \Delta Z_1 \\
+[OFA(1,1-\varphi_1) - OFB(1)] \Delta Z_2\\
+(0)\Delta Z_3
\end{aligned}$$

#### 找数据

需要每个湖的：

- **precipitation** 降水量
- **runoff from the lake's basin**
- **evaporation from the lake surface**
- ice retardation
- diversion
- consumptive use
- groundwater contribution
- **lake surface area**

加粗的是优先寻找的。

waterlevel [waterlevel ](https://www.lre.usace.army.mil/Missions/Great-Lakes-Information/Great-Lakes-Information-2/Water-Level-Data/)

evaporation、precipitation、runoff https://www.glerl.noaa.gov/pubs/tech_reports/glerl-083/UpdatedFiles/

lake surface area https://web.archive.org/web/20120529233616/http://www.epa.gov/glnpo/physfacts.html

consumptive use https://waterusedata.glc.org/data_about_gen.php https://waterusedata.glc.org/graph.php?year=2012&type=summary

#### 代码


$$C_1 = N(1) + QO_{0,m} - OF(1) + QR_1$$

$$c_{1,1} = \dfrac{A_1}{\Delta t} + OFA(1,\varphi_1) + OFB(1) $$
$$c_{1,2} = OFA(1,1-\varphi_1) - OFB(1)$$
$$c_{1,3} = 0$$

注意到数据计算时使用的都是 $m^3$，注意换算单位。


https://docslib.org/doc/12051918/coordinated-great-lakes-regulation-and-routing-model-draft-users

在 user manual 里找到的邪恶参数。

Midlakes 模型大失败。

### 自行 YY

还是使用离散分步累计的方式。

考虑一个湖，由于 elevation 变化极小，可以考虑为表面积不变，那么：

$$\mathrm{d} V = S \times \mathrm{d} Z$$

只需要计算 $\mathrm{d} V$ 也即体积增量。

$$\mathrm{d} V = NBS + flow_{in} - flow_{out}$$

其中：

$$NBS = P - EV + R$$

暂时先考虑了 precipitation 降水，evaporation 蒸发量，runoff to lake 径流量的数据。

为了与实际数据贴合，我们添加一个拟合参量 EX，用于排除各种固定因素的影响：

$$\mathrm{d} V = NBS + flow_{in} - flow_{out} + EX$$

在一个 $\mathrm{d} t$ 的时间增量内进行计算。

那么关键问题就在于，如何计算 $flow_{in}$ 与 $flow_{out}$.

我们依然将湖从上游到下游标号为 $0-4$，那么实际上 $flow_{in}(i) = flow_{out}(i-1)$.

我们实际上只需要 $flow_{out}(0,1,2,3)$ 即可。

其中 $flow_{out}(0)$ 是 St.Mary River 对应的流量，是我们可控制的自变量。

$flow_{out}(4)$ 是最终流入大西洋那条河的流量。

经过简单的式子转化，可以写出：

$$
\begin{aligned}
\mathrm{d} Z_0 &= \dfrac{NBS(0) - flow_{out}(0) + EX(0)}{S_0} \\
\mathrm{d} Z_1 &= \dfrac{NBS(1) + flow_{out}(0) - flow_{out}(1) + EX(1)}{S_1} \\
\mathrm{d} Z_2 &= \dfrac{NBS(2) + flow_{out}(1) - flow_{out}(2) + EX(2)}{S_2} \\
\mathrm{d} Z_3 &= \dfrac{NBS(3) + flow_{out}(2) - flow_{out}(3) + EX(3)}{S_3} \\
\mathrm{d} Z_4 &= \dfrac{NBS(4) + flow_{out}(3) - flow_{out}(4) + EX(4)}{S_4}
\end{aligned}
$$

其中 $flow_{out}(0), flow_{out}(4)$ 可看做已知量。(我们设定的控制量)

则 $\mathrm{d} Z_0$ 可以直接计算得出。

只需在已有的体系中建立 $flow_{out}(1,2,3)$ 的估值函数，即可对该式子进行求解。

事实上，midlakes 模型就是建立了这样的一个估值函数。但是很不幸，估值并不是很成功(

我们认为，流量与如下因素有关：

- 某种线性的相关系数
- 两侧的湖泊的高度差(上面已经有过一定的说明了)
- 两侧的湖泊的加权高度和(借鉴 midlakes)
- 时间(每个月单独分类)
- 某种未知的补偿因子

写出一个毫无根据的式子：

$$flow_{out}(i) = K_{i} \times \left[\varphi_i Z_i - (1-\varphi_i)Z_{i+1} \right]^{a_i} + G_{i} \times \left[ \gamma_i Z_i + (1-\gamma_i) Z_{i+1}\right]^{b_i} + EA_{i}$$

就是干！硬套上去计算(

所有的优化参数：

- K[3]
- varphi[3]
- a[3]
- G[3]
- gamma[3]
- b[3]
- EA[3]
- EX[5]
- STEP

经过了大炼钢铁炼炼炼，出路的最好一路丹也仍然有 $1.4/5$ 的单点数据标准差，重复上个月数据的测试对照组标准差在 $1.68/5$ 级别，以微弱的优势胜出。。。我怀疑是过拟合。

### 未曾设想的道路

使用 Lasso 进行多项式回归，将降水量(P)、蒸发量(EV)、湖面高度(Z)、runoff 量(R)、QO0、QO4(我们可以控制的那两条河的流量)的 $k=[1,2]$ 次方多项式作为自变量，下个月的水位作为因变量，进行多项式回归。

多项式的系数选取是有讲究的。过高的阶数可以在训练数据上更好地拟合，但会出现过拟合的问题。过低的阶数则无法拟合数据。

选取 2008 年至 2020 年中有数据的月份，共 128 个数据，选取后 20 个作为测试集，前 88 个作为训练集。

标准差的计算方法：以当月数据点，预测下个月的五大湖水位，假设当前为第 $t$ 个时间，第 $i$ 个湖预测值为 $p_i$，原值为 $o_i$，则标准差定义为：

$$\sqrt{\dfrac{\sum \limits _t \left(\sum \limits _{i=0}^4 (p_{i,t} - o_{i,t})^2\right) \div 5}{len_T}}$$

用于反映预测下个月数据，对每个湖来说，平均的水位误差情况。

训练集、测试集、全集上的标准差请见后文。

下面是 $k=[1]$ 的拟合情况：


下面是 $k=[1,2]$ 的拟合情况：


下面是 $k=[1,2,3]$ 的拟合情况：


下面是 $k=[1,2,3,4]$ 的拟合情况：


下面是 $k=[1,2,3,4,5]$ 的拟合情况：


下面是 $k$ 的上界不同，对应的训练集、测试集、全集上的标准差：

注意 $k=0$ 相等于取均值(

误差的算术平均值：

(0, 0.31552916441998297, 0.9915595739163559, 0.4121049372051795), 
(1, 0.13319092166831112, 0.27502081091416664, 0.15345233441771905), 
(2, 0.05426636023808525, 0.07052841982003705, 0.056589511606935504), 
(3, 0.04323964290865436, 0.04951125708333195, 0.04413558779075116), 
(4, 0.0423475111699788, 0.04650126258753102, 0.04294090422962909), 
(5, 0.04164210969364411, 0.044440636420801524, 0.0420418992260952)

更大的多项式会导致溢出。结合下文中多步预测的效果，我们最终选择 $k=[1,2,3,4]$ 作为我们的模型使用。

最终选用的模型效果如下：


长期预测检验：我们此前的模型都是单步检验，接下来我们进行长期检验，连续预测 12 个月，使用上一步得到的水位作为下一步的输入。

这个在实际生活中没有什么意义，因为我们肯定是拿最新数据来预测下一步的水位。但是可以检验我们模型的长期稳定性，也可以证明我们的模型用于长期预测也有一定的效力。

time (2019, 0)
delta [-0.0037009351191557016, -0.01197607112308674, 0.02058521252098444, 0.0035795907337501376, -0.05869149920010841]
time (2019, 1)
delta [-0.0327409134708887, -0.005473445800902255, -0.020511316948557123, 0.06403161232515231, 0.06885694172905232]
time (2019, 2)
delta [-0.08911104323007635, -0.020738881641733542, -0.06731304751073708, -0.04608763618952594, -0.01651417669569355]
time (2019, 3)
delta [-0.1359916085512225, -0.0842326587326454, -0.10784900478236636, -0.08069380698356099, -0.17729194298965467]
time (2019, 4)
delta [-0.126862068206691, -0.0958450153925412, -0.10350615783590911, -0.0776411405387023, -0.14385019887636474]
time (2019, 5)
delta [-0.11500831099709785, -0.1027061990506013, -0.09114343940728986, -0.04949046393423373, -0.1266737425602571]
time (2019, 6)
delta [-0.09637167527324664, -0.06415154147560997, -0.060471604154088254, -0.05386418386351011, -0.07153988722529903]
time (2019, 7)
delta [-0.10954518667338675, -0.06569021931974817, -0.04885256816257311, -0.023556777827991482, -0.03500505713671487]
time (2019, 8)
delta [-0.13116355544775615, -0.11479545834092164, -0.015588255019935104, 0.01829884145087135, -0.02329428617225915]
time (2019, 9)
delta [-0.052112952813871516, -0.08842143099764144, 0.023364560997436, 0.05514764350527912, -0.07893285478165524]
time (2019, 10)
delta [-0.05360476348246834, -0.07456811857241519, 0.02514150932856296, 0.044782845054385234, -0.024609391331921415]
time (2019, 11)
delta [-0.06623239432443029, -0.08726533190201735, -0.16522190088610955, -0.1472027387450794, 0.009630621392290095]
time (2020, 1)
delta [-0.017914396445434022, -0.06256256028805751, -0.12261919643950137, -0.13603536041972575, -0.09488741106606824]
time (2020, 2)
delta [-0.036633994806066994, -0.08562330698012488, -0.0977984284302238, -0.10195494618977818, -0.125547728378109]
time (2020, 3)
delta [0.012779721064674732, -0.07453566988738203, -0.060083079820373086, -0.026905659096712498, -0.018850325618572583]
time (2020, 4)
delta [-0.005780227797515636, -0.05875732746085305, -0.020318985196553285, 0.011523548487673452, 0.057340927482940174]
time (2020, 5)
delta [-0.028838887207427888, -0.06198506073945964, -0.031106796343721044, -0.014740972664355922, 0.09079939548657023]
time (2020, 6)
delta [-0.047938580016364085, -0.06302875247288853, -0.03197767637161064, 0.015493560327797695, 0.11050550007369964]
time (2020, 7)
delta [-0.011261000847042624, -0.049713724572853835, -0.04305073212506727, 0.005499195269948132, 0.1596807876277495]
0.07684861290156994

画出预测水位与实际水位的图：


误差的算术平均值在 0.06m 量级。

特别地，这后 20 个月也恰好是我们的测试集。此时模型使用全部数据训练，并进行连续预测。如果只使用训练集，使用连续预测的方式补全测试集，误差就相对较大了，但可以看出走势依然一致：


此时标准差在 0.07m 量级。

是否使用测试集训练，进行连续预测时，误差表现接近。这说明我们的模型并没有发生过拟合。

至此，我们的模型经过了检验，已经证明了其可用性。接下来，我们将基于该五大湖水位模型，建立控制模型。

## 回到第二问

我们的模型已经建立完毕了，接下来是建立控制模型。

控制模型很简单，将 QO0、QO4 作为自变量即可。

此时我们来检验我们的湖水模型关于 QO0、QO4 的灵敏性，避免毫无意义的优化(

取 2019 年 1 月数据，预测 2019 年 2 月情况，改变 QO0、QO4，获得水位变动。

我们控制改变单变量，将 QO0 设置为 $[0, 3 \times QO0]$，用标准差衡量模型输出变化：

$$STD = \sqrt{\dfrac{\sum \limits _{i=0}^4 (now_i - old_i)^2}{5}}$$


让我们将目光聚焦于可以产生调控的那些地方，50% 至 150%：


我们的模型成功将这两个因素纳入考量。

接下来，以 QO0、QO4 的调控倍率作为我们的优化对象，设为 $\alpha$、$\beta$，假设原始流量为 $QO0,QO4$，经过我们调控后的流量为：

$$\alpha \cdot QO0, \beta \cdot QO4$$

在此基础上，结合第一问的目标函数，进行双变量的单目标优化问题求解。

还是按照我们的模型思路，分为单步与多步。

首先看简单的单步调控效果如何吧：自变量为当月的调控因子 $\alpha$、$\beta$，因变量为下个月的预测值，目标最小化下个月的预测值与最佳水位的标准差。

以 2019 年 1 月为例，得出的优化解为：

x0 [1, 1] 0.11868501527225553
Optimization terminated successfully.
         Current function value: 0.105145
         Iterations: 7
         Function evaluations: 24
         Gradient evaluations: 8
[0.57517    0.81219112]

相比于原始值，有约 10% 的提升。

当然，这是一个小数据点。我们对所有的数据进行优化，看看平均优化程度是多少。

算法是：原本对下个月的预测距离最佳水位的标准差为 $original$，改善后为 $now$，则优化率为 $(original-now) \div original$.


当然，还有经典的两曲线对比：


以及改善后的预测值与实际值的对比：


改善的百分比：


可以看到，在绝大部分时间，都可以获得更优的结果。平均而言，可以获得 13.89% 的改善。

单步控制的效果是有限的，因为我们只能控制当月的水坝流量，影响下个月，而一个月的流量并不能大幅改变水位。

接下来，我们进行长期的调控。我们进行为期一年的调控规划，以接下来一年的 $\alpha$、$\beta$ 作为自变量，连续预测一年并计算其与最佳水位的标准差的和。目标函数即调控后的一年预测值与最佳水位的标准差的和。

这是一个 22 个变量的多自变量单目标规划模型，我决定采用遗传算法来求解它。

我们以第一个月为开始，最后一个月为结束，12 个月，前 11 个月的 $\alpha,\beta$ 用于求解下一个月的水位，而最后一个月的将用于求解明年的第一个月，将不纳入考虑。

尝试了使用遗传算法与 SciPy 中的 Minimize 算法进行求解。都可以得到在约束条件下的最优解。

但是 SciPy 中的 SLSQP 算法可以在几秒内完成求解收敛的工作，而遗传算法需要十几分钟。SLSQP 算法在求解这种数学问题上具有非常强大的优势。


上图就是调控后的水位、原始水位、实际水位与目标曲线的曲线。

距离目标水位的差值：


经过计算：

controlled loss 1.6261763184707514
actual loss 2.993322988836306

其中，LOSS 的计算方式是：

假如 $i$ 月份计算得到水位 $Z_{i,0},Z_{i,1},\cdots,Z_{i,4}$，则：

$$LOSS = \sum \limits _{i=0}^{11} \sqrt{\dfrac{\sum \limits _{j=0}^4 (Z_{i,j} - Best_j)^2}{5}}$$

我们的算法对后几个湖有明显的改善效果，对于第一个湖，可能是因为预测本身并不是很准确，所以控制效果不太明显。

总体而言，根据数据，在我们的指标下，我们的控制算法是更优的！

## 第四问

由于缺乏 winter snowpacks、ice jams 等相关数据，这些指标并没有直接被纳入考量。但是降水量是有的。

依然是拿 2017 年举例子，我们将全年的降雨量增大为原先的 $1.3$ 倍，并观察对水位预测的影响：


降水的敏感度：


## 第五问

在我们的第一问中，我们简单地认为所有的湖泊的权利是平等的。如今，我们对湖泊也进行权重处理。

实际上，我们被要求只考虑 Ontario 的利益。

在第二问的 dam control 中，我们修改评价函数，只考虑 Lake Ontario 的评价。修改后，我们重新运行控制算法：



可以看到，我们的 Lake Ontario 直接被锁定到了最佳水位处。强力的控制！

相应的 alpha 与 beta：

[0.96633213 0.59149068 0.3884962  0.17702738 0.4248684  0.97940763
 0.7884512  0.85956903 1.06049942 0.9786624  0.9752747 ] 
 [0.94247665 1.08481515 1.08076656 1.2484085  1.13727732 1.08410161
 1.15640581 0.80239036 0.99035728 1.07056431 1.03231651]

## 可以改进的点

- 使用更有效的水位预测模型
- 细化评价估值，为每年、每个湖、每个月都建立一个单独的评价函数，以确定随时间变化的最佳水位。
- 细化大坝控制函数，考虑大坝控制流速的限制，比如说蓄水量等。
- 将各种环境因素纳入考虑，例如冰塞、地下水、diversion、consumptive use 等。
- 优化模型的长期预测效果。
- 与 RCM(Regional Climate Model) 和 GCM(Global Climate Model) 结合，提供更加详尽的预报数据。
- 结合周边 stakeholders 的量化数据，计算具体的利益考量、更加可靠的评价体系。

## 模型优点缺点

- 优点
- 与midlakes相比，模型直接从数据角度建模，实现简单。只使用了 Lasso 回归，模型复杂度低。
- 对模型进行鲁棒性和灵敏性测试可知，模型具有良好的泛化性能和鲁棒性能，对外界参数变化不敏感，具有一定的自我调节能力。
- 模型运行速度快，几秒即可完成训练到预测的全流程。
- 从模型在2017年数据上的良好变现可知，模型的长期预测准确性较高，具有较强的可实用性。
- 可以将各种环境因素纳入考量。
- 参数多，灵活，可以与外部模型良好交互。