使用穆穆院士提出的CNOP方法对Lorenz63模型进行最优扰动的求解。
1 Lorenz63模型
Lorenz63模型由美国数学家、气象学家爱德华·诺顿·洛伦茨（1917-2008）1963年在论文《Deterministic Nonperiodic Flow》中首次提出。这个三变量无量纲模型简要地描述了大气中常见的热对流现象。
$\frac {dx}{dt}=\sigma (y-x) $
$\frac{dy}{dt}=x(\rho -z)-y$
$\frac{dz}{dt}=xy-\beta z$
其中$x,y,z$为变量，$x$与对流强度成正比例，$y$与上升和下降的温度差成比例，$z$表示平均温度的垂直梯度的偏离。$\sigma,\rho,\beta$为模型参数，$\sigma$是Prandtl数，反映了粘性和热传导之间的比例；$\rho$是Rayleigh数，表示上下层之间的温度差，反映对流强弱；$\beta$表示所研究对流区域的外形比参数。在Lorenz的文章中，这些参数分别取$\sigma=10,\beta=8/3,\rho=28$。
为什么说Lorenz63模型（很多时候人们简称为Lorenz模型）是描述混沌的模型？因为在该模型中，系统初始状态的微小差别会导致迥然不同的状态发展变化，这种就是典型的混沌现象。也就是记者宣传描述的蝴蝶效应（一只蝴蝶在巴西轻拍翅膀，可以导致一个月后德克萨斯州的一场龙卷风）。当我们将该模型计算的结果，绘制在分别以变量$x,y,x$为坐标轴的空间里，图像形如蝴蝶振翅，这或许就是蝴蝶效应的灵感来源。
Lorenz, Edward N., 1963: Deterministic Nonperiodic Flow. Journal of the Atmospheric Sciences, 20, 130-141doi:0.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2020<0130:DNF>2.0.CO;2)
lorenz模式的切线性模式有效性的研究 - 豆丁网
如何通俗地解释混沌理论（chaos）和分岔理论（bifurcation）？
2 模型的差分求解
使用最简单的前向差分：$\frac{dx}{dt}=\frac{x_{i+1}-x_i}{\Delta t}=F(x_i,y_i,z_i)$
于是：$x_{i+1}=x_i+\Delta t F(x_i,y_i,z_i)=x_i+\Delta t \sigma(y_i-x_i)$
其余两个变量的差分类似。
根据前向差分公式，定义一个计算通用的差分计算函数：
1
2
3
| def forward(x,func,h):
x1 = x+h*func(x)
return x1
|
再根据模型右边的具体表达式，定义函数：
1
2
3
4
5
6
7
8
9
| def L63_rhs(x):
import numpy as np
dx=np.ones_like(x);
sigma=10.0; rho=28.0;beta=8/3;
dx[0]=sigma*(x[1]-x[0])
dx[1]=rho*x[0]-x[1]-x[0]*x[2];
dx[2]=x[0]*x[1]-beta*x[2]
return dx
|
以上函数结合构成求解的函数：
1
2
3
4
| def L63_adv_1step(x0,delta_t):
x1=forward(x0,L63_rhs,delta_t)
return x1
|
对函数进行反复迭代调用求解问题：
1
2
3
4
5
6
| time = 1200
delta_t = 0.01
xbase0 = np.array([0,1,0])
xbase = np.zeros([time,3]);xbase[0]=xbase0
for j in range(1,time):
xbase[j] = L63_adv_1step(xbase[j-1],delta_t)
|
绘制三维结果图：
1
2
3
4
5
6
7
8
9
10
11
12
| import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.add_axes(Axes3D(fig))
ax.plot(xbase[:,0], xbase[:,1], xbase[:,2],'r', label='control')
ax.legend()
plt.xlabel('x');plt.ylabel('y');
ax.set_zlabel('z')
|
3 切线性模式的编写和检验
切线性模式近似描述初始扰动随时间的增长规律。切线性模式是非线性模式的一阶近似。
3.1 构造切线性模式
对初始基态$(x,y,z)$添加扰动$(\delta x,\delta y,\delta z)$，得到如下方程：
$\frac{d(x+\delta x)}{dt}=-\sigma(x+\delta x)+\sigma(y+\delta y)$
$\frac{d(y+\delta y)}{dt}=-(x+\delta x)(z+\delta z)+\rho(x+\delta x)-(y+\delta y)$
$\frac{d(z+\delta z)}{dt}=(x+\delta x)(y+\delta y)-\beta (z+\delta z)$
左右约去基态满足的部分（原始Lorenz方程），再省略高阶项（如$\delta x\delta z$等），就得到Lorenz的切线性模式：
$\frac{d\delta x}{dt}=-\sigma \delta x+\sigma \delta y$
$\frac{d\delta y}{dt}=(\rho-z)\delta x-\delta y-x\delta z$
$\frac{d \delta z}{dt}=y\delta x+x\delta y-\beta \delta z$
根据公式写出切线性模式的计算函数：
1
2
3
4
5
6
7
8
9
| def TL_model(x,dx):
dx1 = np.ones_like(dx)
sigma=10.0; rho=28.0;beta=8/3
dx1[0] = sigma*(dx[1]-dx[0])
dx1[1] = (rho-x[2])*dx[0]-dx[1]-x[0]*dx[2]
dx1[2] = x[1]*dx[0]+x[0]*dx[1]-beta*dx[2]
return dx1
|
和原始模式一样，切线性模式也使用前向差分进行迭代计算，注意计算每一步的切线性中，都需要给出相应时刻的基态$x0=[x_i,y_i,z_i]$，这需要我们先正向计算一遍原始模式积分，然后在相应时刻将基本态输入给切线性模式作为”参数“：
1
2
3
4
5
6
7
| def cal_TL(x0,dx0,delta_t):
TL_model2 = partial(TL_model,x0)
dx1=forward(dx0,TL_model2,delta_t)
return dx1
|
3.2 检验切线性模式
切线性模式构造地是否正确，需要进行检验，根据以下公式计算R：
$R=\frac{||\pmb {M}(X_0+\alpha \cdot x_0 )-\pmb {M}(X_0)||}{\alpha ||\rm{\pmb{M}}(X_0) x_0||}=1+O(\alpha)$
其中，$\pmb M$表示非线性模式；$\rm\pmb M$表示$\boldsymbol M$的切线性模式；$X_0$表示参考态的初始场；$x_o$表示$X_0$的初始扰动；$||X||=\sqrt {\sum_ix_i^2}$表示二范数；$0<\alpha<1$。
该检验通过的标准是：
以上两点满足，则说明该切线性模式编写准确。
注意，由于切线性近似只在短距离内满足，所以验证的时候不要迭代非常长的时间，那样效果肯定很差。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
| def check_TL(alpha):
time = 120
x0 = np.array([-6*math.sqrt(2),-6*math.sqrt(2),27])
diff_x0 = x0*0.01
Xtrue = np.zeros([time,3]);Xtrue[0]=x0
Xpert = np.zeros([time,3]);Xpert[0]=x0+alpha*diff_x0
diff_x = np.zeros([time,3]);diff_x[0] = diff_x0
delta_t= 0.01
for j in range(1,time):
Xtrue[j] = L63_adv_1step(Xtrue[j-1],delta_t)
Xpert[j] = L63_adv_1step(Xpert[j-1],delta_t)
diff_x[j] = cal_TL(Xtrue[j-1], diff_x[j-1], delta_t)
print(diff_x[-1])
fmu = alpha*np.linalg.norm(diff_x[-1],ord=2)
print(fmu)
print(abs(Xpert[-1]-Xtrue[-1]))
fzi = np.linalg.norm([Xpert[-1]-Xtrue[-1]],ord=2)
print(fzi)
R = fzi/fmu
print(R)
return R
NN = np.arange(-11,1)
RR = np.ones_like(NN,dtype=float)
for i in range(len(NN)):
n = NN[i]
print(n)
alpha = math.pow(10,n)
RR[i] = check_TL(alpha)
fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(8, 8), dpi=100)
axs.plot(NN, RR, c='k')
axs.set_xlabel("lg(a)", fontdict={'size': 16})
axs.set_ylabel('R', fontdict={'size': 16}, rotation=0)
fig.autofmt_xdate()
plt.show()
|
4 伴随模式的编写和检验
4.1 什么是伴随模式
对于一个内积运算$<,>$（比如最简单的点积），若$\rm{\pmb {M^}}$满足$=<\rm{\pmb {M^}} b,a>$则$\rm{\pmb {M^*}}$是线性模式$\rm{\pmb {M}}$的伴随模式，其中$a,b$均为矢量。
那么伴随模式是如何与梯度联系起来的？
存在一个非线性映射关系$y=M(x)$，若$y_0=M(x_0),y_1=M(x_1)，x_1=x_0+\delta x$，构造$M$的切线性模式$\rm{\pmb M}$使得$\delta y = M(x_1)-M(x_0)=\rm{\pmb M}\delta x$
定义目标函数$J=||y_1-y_0||^2=||\delta y||^2$
则目标函数的梯度为$\delta J=<2\delta y,\delta y>=2<\rm{\pmb M}\delta x,\rm{\pmb M}\delta x>$
若存在$\rm\pmb M^$使得$\delta J=2<\rm{\pmb M}\delta x,\rm{\pmb M}\delta x>=2<\rm{\pmb M^}\rm{\pmb M}\delta x,\delta x>$，则$\rm{\pmb M^*}$是$\rm{\pmb M}$的伴随模式。
考虑梯度的定义：对$x_0$的任一扰动$\delta x$产生的响应函数的一阶变化$\delta J$等于梯度向量$\nabla_{x_0}J$和扰动向量$\delta x$的内积：$\delta J=<\nabla_{x_0}J,\delta x>$
于是，$\nabla_{x_0}J=2\rm{\pmb M^*}\rm{\pmb M}\delta x$
4.2 构造伴随模式并检验
对于切线性模式，伴随模式就是其转置。
根据上面构造的切线性模式，可以写出伴随模式：
$x^=-\sigma x^+(\rho-z)y^+yz^$
$y^=\sigma x^-y^+xz^$
$z^=-xy^-bz^*$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
| def adj_model(xb,xa):
xa0 = np.ones_like(xa)
sigma=10.0; rho=28.0;beta=8/3
xa0[0] = -sigma*xa[0]+(rho-xb[2])*xa[1]+xb[1]*xa[2]
xa0[1] = sigma*xa[0]-xa[1]+xb[0]*xa[2]
xa0[2] = -xb[0]*xa[1]-beta*xa[2]
return xa0
def cal_adj(xb,xa,delta_t):
adj_model2 = partial(adj_model,xb)
xa0 = forward(xa,adj_model2,delta_t)
return xa0
def check_adj():
time = 120
x0 = np.array([-6*math.sqrt(2),-6*math.sqrt(2),27])
diff_x0 = x0*0.01
xbase = np.zeros([time,3]);xbase[0]=x0
diff_x = np.zeros([time,3]);diff_x[0] = diff_x0
delta_t= 0.01
for j in range(1,time):
xbase[j] = L63_adv_1step(xbase[j-1],delta_t)
diff_x[j] = cal_TL(xbase[j-1], diff_x[j-1], delta_t)
xadj = np.zeros([time,3]);xadj[-1]=diff_x[-1]
for j in range(time-1,0,-1):
xadj[j-1] = cal_adj(xbase[j-1], xadj[j], delta_t)
left = np.dot(diff_x[0],xadj[0])
right = np.dot(diff_x[-1],xadj[-1])
print('left=',left)
print('righ=',right)
left = np.dot(diff_x[0],diff_x[0])
xadj0 = cal_adj(xbase[1], diff_x[0], delta_t)
right = np.dot(diff_x[0],xadj0)
print('left=',left)
print('righ=',right)
|
5 目标函数和梯度
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
| def cal_J(x0,dx0):
time = 120
xbase = np.zeros([time,3]);xbase[0]=x0
xdiff = np.zeros([time,3]);xdiff[0]=x0+dx0
delta_t= 0.01
for j in range(1,time):
xbase[j] = L63_adv_1step(xbase[j-1],delta_t)
xdiff[j] = L63_adv_1step(xdiff[j-1],delta_t)
diff_xn = xdiff[-1]-xbase[-1]
J = -pow(np.linalg.norm(diff_xn,ord=2),2)
return J,diff_xn,xdiff
def cal_G(x0,dx0):
time = 120
delta_t = 0.01
J,diff_xn,xdiff = cal_J(x0,dx0)
xadj = np.zeros([time,3]);xadj[-1]=diff_xn
for j in range(time-1,0,-1):
xadj[j-1] = cal_adj(xdiff[j-1], xadj[j], delta_t)
G = -2*xadj[0]
print('G=',G)
return G
def check_G(alpha,x0,dx0):
delta_t = 0.01
Jbase,diff_x,xdiff = cal_J(x0,dx0)
G = cal_G(x0,dx0)
normG = np.linalg.norm(G,ord=2)
print(normG)
dxalpha = dx0 + alpha*G/normG
Jalpha = cal_J(x0,dxalpha)[0]
print(Jalpha)
R = (Jalpha - Jbase)/(alpha*normG)
return R
|
6 SPG2优化算法的编写
6.1 SPG2算法介绍
6.2 编写SPG2算法
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
| def func_project(x_centre,x_out):
diff_x = x_out - x_centre
radius_lim = 0.01*np.linalg.norm(x_centre,ord=2)
radius_diff = np.linalg.norm(diff_x,ord=2)
if radius_diff <= radius_lim:
x_pro = x_out
else:
x_pro = diff_x/radius_diff*radius_lim+x_centre
return x_pro
time = 1200
delta_t = 0.01
xbase0 = np.array([-6*math.sqrt(2),-6*math.sqrt(2),27])
xdiff0 = xbase0+np.random.rand(3,)
xdiff0 = func_project(xbase0, xdiff0)
diff_x0 = xdiff0-xbase0
maxit = 1000
maxfun = 2000
k = 0
gamma = 1e-4
alpha_spg0 = 1.0
alpha_max = 1e3
alpha_min = 1e-3
sigma1 = 0.1
sigma2 = 0.9
g0 = cal_G(xbase0,diff_x0)
d0 = func_project(xbase0, xdiff0-alpha_spg0*g0)-xdiff0
J = []
J0 = cal_J(xbase0,diff_x0)[0]
J.append(J0)
JP = []
def spg_end_check(x0,x,g):
eps = 0
xp = abs(func_project(x0,x-g)-x)
if max(xp) <= eps:
print('SPG2 end!,info = 0')
return True
else:
return False
m = 1
while spg_end_check(xbase0,xdiff0,g0) == False:
m = m+1
if m >= maxit:
break
else:
print('m=',m)
lmd = 1
xplus = xdiff0+d0
delta = np.dot(g0,d0)
Jplus = cal_J(xbase0,xplus-xbase0)[0]
JP.append(Jplus)
while (Jplus > (max(max(J),Jplus)+gamma*lmd*delta)):
lmd_temp = -0.5*(lmd**2)*delta/(Jplus-J0-lmd*delta)
if lmd_temp >= sigma1 and lmd_temp <= (sigma2*lmd):
lmd = lmd_temp
break
else:
lmd = 0.5*lmd
xplus = xdiff0+lmd*d0
Jplus = cal_J(xbase0,xplus-xbase0)[0]
lmd0 = lmd
xdiff1 = xdiff0+lmd0*d0
g1 = cal_G(xbase0,xdiff1-xbase0)
s0 = xdiff1-xdiff0
y0 = g1-g0
beta0 = np.dot(s0,y0)
if beta0 <= 0:
alpha_spg1 = alpha_max
else:
alpha_spg1 = min(alpha_max,max(alpha_min,np.dot(s0,s0)/beta0))
xdiff0 = xdiff1
alpha_spg0 = alpha_spg1
diff_x0 = xdiff0-xbase0
g0 = cal_G(xbase0,diff_x0)
d0 = func_project(xbase0, xdiff0-alpha_spg0*g0)-xdiff0
J0 = cal_J(xbase0,diff_x0)[0]
J.append(J0)
|