随机微分方程（Stochastic Differential Equations，SDEs）是描述系统在随机扰动下动态行为的数学模型。与常微分方程（Ordinary Differential Equations，ODEs）不同，SDEs考虑了系统中随机性或不确定性，通常用于建模物理、金融、生物等领域的动态系统。以下是对随机微分方程的详细介绍。

1. 随机微分方程的基本形式

随机微分方程通常写作如下形式：

• Xt：表示在时间ttt的状态（通常是一个随机过程）。
• μ(Xt,t)：漂移项，表示系统的确定性部分，描述了状态变量随时间的变化趋势。
• σ(Xt,t)：扩散项，表示随机性或噪声的影响程度。
• Wt：表示标准布朗运动（或维纳过程），是描述随机扰动的基础。

2. 随机微分方程的解释

• 漂移项（Drift Term）：漂移项μ(Xt,t)\mu(X_t, t)μ(Xt,t)控制系统的平均行为。它是系统在没有噪声影响时的演变趋势。
• 扩散项（Diffusion Term）：扩散项σ(Xt,t)\sigma(X_t, t)σ(Xt,t)控制系统在随机扰动下的波动程度。它反映了噪声对系统行为的影响。
• 布朗运动（Brownian Motion）：布朗运动是随机过程的一种，表示连续时间内的随机变化，具有独立增量和连续轨迹的特性。它是SDE中的主要随机成分。

3. 解决随机微分方程

解决SDE的方法与常微分方程不同，常用的方法有：

• 伊藤积分（Itô Integral）：伊藤积分是一种特殊类型的积分，适用于随机过程，常用于求解SDE。
• 伊藤公式（Itô's Lemma）：类似于常微分方程中的链式法则，伊藤公式用于处理SDE的求导和积分，常用于构造解的形式。
• 数值解法：常用的方法包括： 欧拉-马鲁扬方法（Euler-Maruyama Method）：用于求解SDE的显式数值方法。 随机Runge-Kutta方法：扩展了常规Runge-Kutta方法以适应随机情况。

4. 应用领域

随机微分方程在多个领域有广泛应用：

• 金融工程：用于建模股票价格、利率、期权定价等金融市场中的随机现象。著名的Black-Scholes模型就是基于SDE构建的。
• 物理学：用于描述粒子的随机运动，例如布朗运动等现象。
• 生物学：用于建模种群动态、遗传变异等生物现象中的随机过程。
• 控制理论：在控制系统中引入随机扰动的影响，以设计更鲁棒的控制器。

5. 示例

下面是一个简单的随机微分方程示例：

• 这个方程描述了某种量（例如某种资源或价格）如何在平均值μ周围波动，θ是收敛速度，σ是波动率。

6. Julia中的实现

如前所示，可以使用Julia的包来求解SDE。该库提供了强大的功能，可以方便地定义和解决随机微分方程。

代码块

using DifferentialEquations

function sde_system!(du, u, p, t)
   du[1] = ...
   # 其他的方程
end

# 设置问题
sde_prob = SDEProblem(sde_system!, u0, tspan, p)
sol = solve(sde_prob)

      复制成功
     
     
     
     

在科学计算领域，特别是ODE类微分方程的求解，Julia已经实现并覆盖了最大部分的求解算法，相比于其他科学计算软件（MATLAB）

图片来源：https://www.stochasticlifestyle.com/comparison-differential-equation-solver-suites-matlab-r-julia-python-c-fortran/#:~:text=For%20the%20current%20state%20of%20the%20reproducible%20benchmarks

视频里面up是用matlab的ode45求解器实现的，由于原论文的控制器参数不能正常复现实验结果，所以参考了视频中的参数设置，完整julia代码如下

using DifferentialEquations
using Plots

function sig(x, alpha)
    return abs(x)^alpha * sign(x)
end

function d_t(t)
    return 2 - cos(t)
end

function control_law(x1, x2, x3)
    return -35.9 * x1 - 48.15 * x2 - 35.47 * x3
end

function sde01!(dx, x, p, t)
    u = control_law(x[1], x[2], x[3])
    dx[1] = d_t(t) * sig(x[2], 5/7) + 1/10*sin(x[3])
    dx[2] = sig(x[3], 5/7)
    dx[3] = sig(u, 5/7)
end

function σ_sde01!(dx, x, p, t)
    dx[1] = 1/6*sig(x[3],6/7)
    dx[2] = 0.0
    dx[3] = 0.0
end

prob = SDEProblem(sde01!, σ_sde01!, [0.5, -2.0, 3.0], [0.0, 6.0])
sol = solve(prob)
plot(sol, idx=(1, 2, 3))

function control_law(x)
    x1, x2, x3 = x
    return -35.9 * x1 - 48.15 * x2 - 35.47 * x3
end

plot(sol.t, control_law.(sol.u))

      复制成功
     
     
     
     

states of x

u

与视频中的求解结果一致。

注：Julia支持重载，当参数列表不同时为不同的函数，调用时根据实际情况而定

7. Python中的实现

同时，下面给出python的求解，也适用于其它ode方程的初值问题，

代码块

import numpy as np 
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 定义控制律 u
def control_law(x1, x2, x3):
    return -35.9 * x1 - 48.15 * x2 - 35.47 * x3

def sig(x, alpha):
    return np.abs(x)**alpha*np.sign(x)

# 定义随机微分方程 (SDE)
def sde_system(t, x):
    x1, x2, x3 = x
    
    # 定义 d(t)
    d_t = 2 - np.cos(t)
    
    # 控制律 u
    u = control_law(x1, x2, x3)
     # 随机扰动部分
    noise =  np.random.normal()
    # 微分方程的确定性部分
    dx1_dt = d_t * sig(x2, 5/7) + (1/10) * np.sin(x3) + (1/6) * sig(x3,6/7) * noise
    dx2_dt = sig(x3, 5/7)
    dx3_dt = sig(u, 5/7)
    
   
    
    return np.array([dx1_dt, dx2_dt, dx3_dt])

# 设置初始条件和时间范围
x0 = [0.5, -2, 3]  # 初始条件 [x1(0), x2(0), x3(0)]
T = 10   # 时间范围
t_span = [0, T]
t_eval = np.linspace(0, 10, 5000)  # 用于评估的时间点

# 噪声强度
noise_strength = 0.1

# 使用 solve_ivp 来求解 SDE
sol = solve_ivp(sde_system, t_span, x0, t_eval=t_eval, method='RK45')

# 提取解
x1_sol = sol.y[0]
x2_sol = sol.y[1]
x3_sol = sol.y[2]

# 计算控制律 u(t) 在每个时间点的值
u_sol = control_law(x1_sol, x2_sol, x3_sol)

# 绘制状态变量 x1, x2, x3 和控制律 u 的结果
plt.figure(figsize=(12, 8))

# 状态变量图
plt.subplot(2, 1, 1)
plt.plot(t_eval, x1_sol, label="x1(t)", color="r")
plt.plot(t_eval, x2_sol, label="x2(t)", color="g")
plt.plot(t_eval, x3_sol, label="x3(t)", color="b")
plt.xlim([0,5])
plt.xlabel('Time (t)')
plt.ylabel('States (x1, x2, x3)')
plt.title('Closed-loop System Response')
plt.legend()
plt.grid(True)

# 控制律 u 图
plt.subplot(2, 1, 2)
plt.plot(t_eval, u_sol, label="u(t)", color="m")
plt.xlim([0,6])
plt.xlabel('Time (t)')
plt.ylabel('Control Law u(t)')
plt.title('Control Law u(t) over Time')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()



      复制成功
     
     
     
     

总结

随机微分方程是理解和建模具有随机性的动态系统的重要工具。通过结合确定性和随机性，它们能够提供对现实世界中许多现象的深入理解和预测能力。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删