当前位置：首页 > article >正文

DNN(深度神经网络)近似 Lyapunov 函数

article 2025/2/7 15:44:09

import torch
import torch.nn as nn
import torch.optim as optim
import matplotlib.pyplot as plt
# from torchviz import make_dot
import torchviz


# 1. Lyapunov 函数近似器（MLP 结构）
class LyapunovNet(nn.Module):
    def __init__(self, input_dim, hidden_dim=32):
        super(LyapunovNet, self).__init__()
        self.model = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, 1)  # 输出一个标量，表示Lyapunov函数值
        )
    
    def forward(self, x):
        return self.model(x)

# 2. 梯度伴随网络（计算 Lyapunov 函数的梯度）
def compute_gradient(model, x):
    x.requires_grad_(True)
    V = model(x)
    grad_V = torch.autograd.grad(V.sum(), x, create_graph=True)[0]
    return V, grad_V

# 3. 训练数据（随机生成一些状态数据）
def generate_data(num_samples=1000, state_dim=2):
    return torch.randn(num_samples, state_dim) * 5  # 扩大范围

# 4. 网络权重初始化
def init_weights(m):
    if isinstance(m, nn.Linear):
        torch.nn.init.xavier_uniform_(m.weight)
        if m.bias is not None:
            torch.nn.init.zeros_(m.bias)

# 5. 训练 Lyapunov 网络并绘制损失曲线
def train_lyapunov_net(state_dim=2, epochs=500, lr=0.001):
    model = LyapunovNet(input_dim=state_dim)
    model.apply(init_weights)  # 重新初始化网络
    optimizer = optim.Adam(model.parameters(), lr=lr)
    loss_history = []
    
    for epoch in range(epochs):
        x = generate_data()
        V, grad_V = compute_gradient(model, x)
        
        # 修正损失函数
        loss = torch.mean(torch.relu(-grad_V.sum(dim=1))) + torch.mean(torch.relu(-V))
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        
        loss_history.append(loss.item())
        
        if epoch % 50 == 0:
            print(f"Epoch {epoch}, Loss: {loss.item():.6f}")
    
    # 绘制损失曲线
    plt.plot(loss_history)
    plt.xlabel("Epoch")
    plt.ylabel("Loss")
    plt.title("Training Loss Curve")
    plt.grid()
    plt.show()
    
    return model

# 6. 可视化神经网络计算图
def visualize_lyapunov_network(model, state_dim=2):
    x = torch.randn(1, state_dim, requires_grad=True)  # 生成一个测试输入
    V = model(x)
    dot = make_dot(V, params=dict(model.named_parameters()))
    dot.format = 'png'
    dot.render('lyapunov_network')  # 生成 PNG 图片
    dot.view()  # 打开图像

if __name__ == "__main__":
    trained_model = train_lyapunov_net()
    visualize_lyapunov_network(trained_model)

在这里插入图片描述

这段代码的核心目标是使用神经网络近似 Lyapunov 函数，并通过梯度信息优化其参数，以确保 Lyapunov 函数在物理系统的状态空间中满足稳定性条件。以下是详细的解析：

1. 代码整体架构

该代码主要包含 6 个部分：

定义 Lyapunov 函数近似器（MLP 神经网络）。
计算 Lyapunov 函数的梯度（用于优化）。
生成训练数据（模拟状态空间点）。
初始化网络权重（Xavier 初始化）。
训练 Lyapunov 网络（基于梯度信息优化 Lyapunov 函数）。
可视化网络结构（绘制计算图）。

2. 详细解析每个部分

(1) Lyapunov 函数近似器

class LyapunovNet(nn.Module):
    def __init__(self, input_dim, hidden_dim=32):
        super(LyapunovNet, self).__init__()
        self.model = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, 1)  # 输出一个标量，表示Lyapunov函数值
        )
    
    def forward(self, x):
        return self.model(x)

📌 作用：

这是一个

多层感知机（MLP）

结构：
- 输入 x 代表系统状态变量（维度 input_dim）。
- 两个隐藏层，每层 hidden_dim=32，激活函数为 ReLU。
- 输出是一个标量，表示 Lyapunov 函数值 V(x)。

📌 数学意义：

训练后，神经网络会学习到一个 Lyapunov 函数近似 V(x)，用于判定系统稳定性。

(2) 计算 Lyapunov 函数的梯度

def compute_gradient(model, x):
    x.requires_grad_(True)  # 使 x 可微分
    V = model(x)  # 计算 Lyapunov 函数值
    grad_V = torch.autograd.grad(V.sum(), x, create_graph=True)[0]  # 计算梯度
    return V, grad_V

📌 作用：

计算 V(x) 对 x 的梯度 ∇V(x)。
torch.autograd.grad() 自动计算梯度，用于 Lyapunov 函数优化。

📌 数学意义：

$V (x)$ 的梯度 $\nabla V (x)$

需要满足：

$\frac{dV}{dt} = \nabla V(x) \cdot f(x) \leq0$

这个条件保证 Lyapunov 函数单调递减，从而保证系统稳定性。

(3) 生成训练数据

def generate_data(num_samples=1000, state_dim=2):
    return torch.randn(num_samples, state_dim) * 5  # 生成随机状态点

📌 作用：

生成服从标准正态分布的随机状态数据，用于训练 Lyapunov 神经网络。

📌 数学意义：

这些数据点可以看作是随机抽取的系统状态，用于训练神经网络使其学会合适的 Lyapunov 函数。

(4) 网络权重初始化

辑def init_weights(m):
    if isinstance(m, nn.Linear):
        torch.nn.init.xavier_uniform_(m.weight)
        if m.bias is not None:
            torch.nn.init.zeros_(m.bias)

📌 作用：

采用 Xavier（Glorot）初始化 线性层权重，保证网络初始梯度不会过大或过小。
bias 初始化为 0。

📌 数学意义：

Xavier 初始化

$\sim U\left(-\sqrt{\frac{6}{n_{in} + n_{out}}}, \sqrt{\frac{6}{n_{in} + n_{out}}}\right)$
- 其中 $n_{in}$ 和 $n_{out}$ 分别是输入和输出的神经元数目。

(5) 训练 Lyapunov 网络

def train_lyapunov_net(state_dim=2, epochs=500, lr=0.001):
    model = LyapunovNet(input_dim=state_dim)
    model.apply(init_weights)  # 重新初始化网络
    optimizer = optim.Adam(model.parameters(), lr=lr)
    loss_history = []
    
    for epoch in range(epochs):
        x = generate_data()
        V, grad_V = compute_gradient(model, x)
        
        # 修正损失函数
        loss = torch.mean(torch.relu(-grad_V.sum(dim=1))) + torch.mean(torch.relu(-V))
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        
        loss_history.append(loss.item())
        
        if epoch % 50 == 0:
            print(f"Epoch {epoch}, Loss: {loss.item():.6f}")
    
    # 绘制损失曲线
    plt.plot(loss_history)
    plt.xlabel("Epoch")
    plt.ylabel("Loss")
    plt.title("Training Loss Curve")
    plt.grid()
    plt.show()
    
    return model

📌 作用：

优化 Lyapunov 函数 V(x)，使其满足 Lyapunov 稳定性条件

：
- V(x) > 0（Lyapunov 函数应该是正的）。
- ∇V(x) ≤ 0（Lyapunov 函数应该单调递减）。

📌 数学意义：

损失函数

$L = E [ma x (0, - \nabla V (x))] + E [ma x (0, - V (x))]$
- 其中：
  - torch.relu(-grad_V.sum(dim=1))：强制 ∇V(x) ≤ 0，确保 Lyapunov 单调递减。
  - torch.relu(-V)：确保 V(x) ≥ 0。

(6) 可视化 Lyapunov 神经网络

def visualize_lyapunov_network(model, state_dim=2):
    x = torch.randn(1, state_dim, requires_grad=True)  # 生成一个测试输入
    V = model(x)
    dot = make_dot(V, params=dict(model.named_parameters()))
    dot.format = 'png'
    dot.render('lyapunov_network')  # 生成 PNG 图片
    dot.view()  # 打开图像

📌 作用：

使用 torchviz 绘制计算图，直观展示神经网络结构。

3. main-run

if __name__ == "__main__":
    trained_model = train_lyapunov_net()
    visualize_lyapunov_network(trained_model)

📌 作用：

训练 Lyapunov 网络。
可视化网络结构。

4. 总结

模块	作用
`LyapunovNet`	使用 MLP 近似 `V(x)`
`compute_gradient`	计算 `∇V(x)`
`generate_data`	生成随机状态数据
`init_weights`	Xavier 权重初始化
`train_lyapunov_net`	训练神经网络，使其满足 Lyapunov 条件
`visualize_lyapunov_network`	可视化计算图

这个代码的核心思想是通过深度学习构建 Lyapunov 函数，确保 V(x) 单调递减，从而评估系统的瞬态稳定性。🔥

5.延伸

将以上的损失函数改为平方-以下为改后的损失函数公式展示

给定状态 $x$ 的动态系统：

$\frac{dx}{dt} = f(x)$

定义 Lyapunov 近似函数 $V (x)$ 及其梯度：

$\nabla V = \frac{\partial V}{\partial x}$

损失函数由两个部分组成：

Lyapunov 下降约束（确保 $\frac{dV}{dt} \leq 0$ ）：

$\mathcal{L}_{\text{decay}} = \mathbb{E}_x \left[ \max\left( 0, -\nabla V \cdot f(x) \right) \right]$

其中 $\nabla V \cdot f(x) = \sum_{i} \frac{\partial V}{\partial x_i} f_i(x)$ 是 Lyapunov 函数关于时间的导数。
Lyapunov 函数非负性约束（确保 $V (x) \geq 0$ ）：

$\mathcal{L}_{\text{positivity}} = \mathbb{E}_x \left[ V(x)^2 \right]]$

这里我们使用平方损失 $V(x)^2$ 代替 ReLU 约束 $ma x (0, - V (x))$ 以增强稳定性。

总损失函数：

$\mathcal{L}_{\text{decay}} + \mathcal{L}_{\text{positivity}}$

该损失函数的优化目标是找到一个满足 Lyapunov 条件的近似函数 $V (x)$ ，使得：

Lyapunov 下降条件： $\frac{dV}{dt} \leq 0$ 适用于所有 $x$ 。
Lyapunov 函数非负性： $V (x) \geq 0$ 确保稳定性分析的正确性。

这样，我们训练出的神经网络可以逼近一个合适的 Lyapunov 函数，从而用于评估瞬态稳定性

代码

import torch
import torch.nn as nn
import torch.optim as optim
import matplotlib.pyplot as plt
from torchviz import make_dot

# 1. Lyapunov 函数近似器（MLP 结构）
class LyapunovNet(nn.Module):
    def __init__(self, input_dim, hidden_dim=32):
        super(LyapunovNet, self).__init__()
        self.model = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, 1)  # 输出一个标量，表示Lyapunov函数值
        )
    
    def forward(self, x):
        return self.model(x)

# 2. 定义系统动力学方程（如果已知）
def system_dynamics(x):
    return -x  # 示例：简单的线性动力学 dx/dt = -x

# 3. 计算 Lyapunov 函数的梯度
def compute_gradient(model, x):
    x.requires_grad_(True)
    V = model(x)
    grad_V = torch.autograd.grad(V.sum(), x, create_graph=True)[0]
    return V, grad_V

# 4. 训练数据（随机生成一些状态数据）
def generate_data(num_samples=1000, state_dim=2):
    return torch.randn(num_samples, state_dim) * 5  # 扩大范围

# 5. 网络权重初始化
def init_weights(m):
    if isinstance(m, nn.Linear):
        torch.nn.init.xavier_uniform_(m.weight)
        if m.bias is not None:
            torch.nn.init.zeros_(m.bias)

# 6. 训练 Lyapunov 网络并绘制损失曲线
def train_lyapunov_net(state_dim=2, epochs=500, lr=0.001):
    model = LyapunovNet(input_dim=state_dim)
    model.apply(init_weights)  # 重新初始化网络
    optimizer = optim.Adam(model.parameters(), lr=lr)
    loss_history = []
    
    for epoch in range(epochs):
        x = generate_data()
        V, grad_V = compute_gradient(model, x)
        dx_dt = system_dynamics(x)  # 获取系统的动态方程
        
        # 1. 约束 Lyapunov 下降条件: dV/dt = ∇V · f(x) ≤ 0
        lyapunov_decay = torch.mean(torch.relu(-(grad_V * dx_dt).sum(dim=1)))
        # 2. 确保 Lyapunov 函数 V(x) 是非负的
        lyapunov_positivity = torch.mean(V**2)  # 使 V(x) 逼近 0 而不是仅仅大于 0
        
        loss = lyapunov_decay + lyapunov_positivity
        
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        
        loss_history.append(loss.item())
        
        if epoch % 50 == 0:
            print(f"Epoch {epoch}, Loss: {loss.item():.6f}")
    
    # 绘制损失曲线
    plt.plot(loss_history)
    plt.xlabel("Epoch")
    plt.ylabel("Loss")
    plt.title("Training Loss Curve")
    plt.grid()
    plt.show()
    
    return model

# 7. 可视化神经网络计算图
def visualize_lyapunov_network(model, state_dim=2):
    x = torch.randn(1, state_dim, requires_grad=True)  # 生成一个测试输入
    V = model(x)
    dot = make_dot(V, params=dict(model.named_parameters()))
    dot.format = 'png'
    dot.render('lyapunov_network')  # 生成 PNG 图片
    dot.view()  # 打开图像

if __name__ == "__main__":
    trained_model = train_lyapunov_net()
    visualize_lyapunov_network(trained_model)