用Python可视化李雅普诺夫函数从数学定义到动态理解在控制理论和非线性系统分析中李雅普诺夫函数就像是一盏明灯为我们判断系统稳定性提供了强有力的工具。但传统的教科书式讲解往往让初学者陷入数学公式的泥潭难以建立直观理解。本文将带你用Python的可视化工具把抽象的V(x)正定、导函数负定转化为生动的动态图形让稳定性分析变得触手可及。1. 李雅普诺夫函数的核心思想可视化李雅普诺夫函数本质上是一个能量函数——它描述系统在状态空间中的能量变化。想象一个小球在碗底滚动碗的形状就是我们的李雅普诺夫函数而小球的运动轨迹展示了系统状态的变化。用Python实现这个比喻非常简单import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 创建碗状曲面 - 这就是我们的Lyapunov函数V(x) x np.linspace(-5, 5, 100) y np.linspace(-5, 5, 100) X, Y np.meshgrid(x, y) V X**2 Y**2 # 一个简单的二次型Lyapunov函数 # 绘制3D图形 fig plt.figure(figsize(10, 7)) ax fig.add_subplot(111, projection3d) ax.plot_surface(X, Y, V, cmapviridis, alpha0.8) # 添加小球运动轨迹模拟 t np.linspace(0, 10, 100) x_traj 4 * np.exp(-0.5 * t) * np.cos(t) y_traj 4 * np.exp(-0.5 * t) * np.sin(t) V_traj x_traj**2 y_traj**2 ax.plot(x_traj, y_traj, V_traj, r-, linewidth2, label系统轨迹) ax.scatter(x_traj[0], y_traj[0], V_traj[0], colorblue, s100, label初始状态) ax.scatter(0, 0, 0, colorgreen, s100, label平衡点) ax.set_xlabel(状态变量x1) ax.set_ylabel(状态变量x2) ax.set_zlabel(V(x)) ax.legend() plt.title(李雅普诺夫函数与系统轨迹) plt.show()这段代码展示了几个关键概念V(x)碗状曲面表示我们的Lyapunov函数红色轨迹系统状态随时间变化从高处逐渐向平衡点(碗底)移动能量衰减轨迹在V(x)上的投影高度不断降低直观展示了dV/dt 0提示运行这段代码时可以尝试旋转3D图形从不同角度观察系统如何滑向平衡点。2. 相平面分析与Lyapunov函数导数理解Lyapunov函数的导数dV/dt同样重要。在相平面中我们可以同时绘制系统轨迹和Lyapunov函数的等高线def system_ode(x, t): 定义一个简单的稳定系统 x1, x2 x dx1dt -x1 x2 dx2dt -x1 - x2 return [dx1dt, dx2dt] from scipy.integrate import odeint # 计算轨迹 t np.linspace(0, 10, 100) initial_conditions [[2, 2], [-3, 1], [1, -2]] trajectories [odeint(system_ode, ic, t) for ic in initial_conditions] # 计算Lyapunov函数及其导数 def V(x1, x2): return x1**2 x2**2 def dVdt(x1, x2): return 2*x1*(-x1 x2) 2*x2*(-x1 - x2) # 绘制相平面 plt.figure(figsize(12, 6)) x1 np.linspace(-4, 4, 20) x2 np.linspace(-4, 4, 20) X1, X2 np.meshgrid(x1, x2) # 绘制Lyapunov函数等高线 contour plt.contour(X1, X2, V(X1, X2), levels10, colorsgray, alpha0.5) plt.clabel(contour, inlineTrue, fontsize8) # 绘制轨迹 colors [r, b, g] for traj, color in zip(trajectories, colors): plt.plot(traj[:, 0], traj[:, 1], color, linewidth2) plt.scatter(traj[0, 0], traj[0, 1], colorcolor, s50, labelf初始点({traj[0, 0]}, {traj[0, 1]})) # 绘制dV/dt的向量场 dVdt_x1 2*X1*(-X1 X2) dVdt_x2 2*X2*(-X1 - X2) plt.quiver(X1, X2, dVdt_x1, dVdt_x2, scale100, colorpurple, alpha0.6, labeldV/dt方向) plt.xlabel(x1) plt.ylabel(x2) plt.title(相平面中的Lyapunov函数分析) plt.legend() plt.grid(True) plt.show()这个可视化展示了几个关键点等高线表示Lyapunov函数V(x)的值越靠近中心(平衡点)值越小轨迹不同初始条件下的系统状态演化都趋向于平衡点向量场展示了dV/dt的方向始终指向V(x)减小的方向3. 交互式Lyapunov函数探索静态图像有时难以完全展示Lyapunov函数的动态特性。使用Plotly可以创建交互式可视化import plotly.graph_objects as go # 创建交互式3D图形 fig go.Figure(data[ go.Surface(zV, xX, yY, colorscaleViridis, opacity0.8), go.Scatter3d(xx_traj, yy_traj, zV_traj, modelines, linedict(colorred, width4), name系统轨迹) ]) # 添加初始点和平衡点 fig.add_trace(go.Scatter3d(x[x_traj[0]], y[y_traj[0]], z[V_traj[0]], modemarkers, markerdict(size8, colorblue), name初始状态)) fig.add_trace(go.Scatter3d(x[0], y[0], z[0], modemarkers, markerdict(size8, colorgreen), name平衡点)) # 设置图形布局 fig.update_layout( title交互式Lyapunov函数可视化, scenedict( xaxis_titlex1, yaxis_titlex2, zaxis_titleV(x), cameradict( eyedict(x1.5, y1.5, z0.8) ) ), width800, height600 ) fig.show()交互式可视化允许你旋转图形从任意角度观察缩放查看细节悬停查看具体数值点击图例显示/隐藏特定元素4. 非线性系统的Lyapunov函数分析前面的例子都是线性系统现在让我们看一个非线性系统的例子——著名的Van der Pol振荡器def van_der_pol(x, t, mu1.0): x1, x2 x dx1dt x2 dx2dt mu*(1 - x1**2)*x2 - x1 return [dx1dt, dx2dt] # 定义Lyapunov函数候选 def V_vdp(x1, x2): return 0.5*(x1**2 x2**2) # 计算轨迹 t np.linspace(0, 20, 1000) initial_conditions [[0.5, 0.5], [2, 0], [-1, -1]] trajectories_vdp [odeint(van_der_pol, ic, t) for ic in initial_conditions] # 绘制相平面和Lyapunov函数分析 plt.figure(figsize(14, 6)) # 相平面 plt.subplot(1, 2, 1) x1 np.linspace(-3, 3, 20) x2 np.linspace(-3, 3, 20) X1, X2 np.meshgrid(x1, x2) # 计算向量场 dx1dt X2 dx2dt 1*(1 - X1**2)*X2 - X1 plt.streamplot(X1, X2, dx1dt, dx2dt, colorgray, density1.5) # 绘制轨迹 colors [r, b, g] for traj, color in zip(trajectories_vdp, colors): plt.plot(traj[:, 0], traj[:, 1], color, linewidth2) plt.scatter(traj[0, 0], traj[0, 1], colorcolor, s50) plt.xlabel(x1) plt.ylabel(x2) plt.title(Van der Pol振荡器相平面) plt.grid(True) # Lyapunov函数导数分析 plt.subplot(1, 2, 2) dVdt_vdp X1*X2 X2*(1*(1 - X1**2)*X2 - X1) contour plt.contourf(X1, X2, dVdt_vdp, levels20, cmapcoolwarm) plt.colorbar(contour, labeldV/dt) plt.contour(X1, X2, V_vdp(X1, X2), levels10, colorsblack, alpha0.3) plt.xlabel(x1) plt.ylabel(x2) plt.title(Lyapunov函数导数(dV/dt)分析) plt.grid(True) plt.tight_layout() plt.show()这个例子展示了非线性系统的复杂相平面极限环的存在Lyapunov函数导数在不同区域的表现右图的热图展示了dV/dt的值分布蓝色区域dV/dt 0 (系统能量在减少)红色区域dV/dt 0 (系统能量在增加)注意对于Van der Pol振荡器我们选择的V(x)并不是一个真正的Lyapunov函数因为dV/dt在某些区域为正。这正好展示了寻找合适Lyapunov函数的挑战性。5. 高级可视化技巧动画展示为了让理解更加直观我们可以创建动画来展示Lyapunov函数随时间的变化from matplotlib.animation import FuncAnimation from IPython.display import HTML # 设置图形 fig, (ax1, ax2) plt.subplots(1, 2, figsize(14, 6)) # 准备轨迹数据 traj trajectories[0] # 使用第一个轨迹 v_values V(traj[:, 0], traj[:, 1]) # 相平面设置 ax1.set_xlim(-4, 4) ax1.set_ylim(-4, 4) ax1.set_xlabel(x1) ax1.set_ylabel(x2) ax1.set_title(相平面轨迹) ax1.grid(True) contour ax1.contour(X1, X2, V(X1, X2), levels10, colorsgray, alpha0.5) line, ax1.plot([], [], r-, linewidth2) point, ax1.plot([], [], ro, markersize8) # Lyapunov函数值设置 ax2.set_xlim(0, len(t)) ax2.set_ylim(0, max(v_values)*1.1) ax2.set_xlabel(时间步) ax2.set_ylabel(V(x)) ax2.set_title(Lyapunov函数值随时间变化) ax2.grid(True) v_line, ax2.plot([], [], b-, linewidth2) def init(): line.set_data([], []) point.set_data([], []) v_line.set_data([], []) return line, point, v_line def animate(i): # 更新相平面 line.set_data(traj[:i, 0], traj[:i, 1]) point.set_data(traj[i, 0], traj[i, 1]) # 更新Lyapunov函数值 v_line.set_data(range(i), v_values[:i]) return line, point, v_line ani FuncAnimation(fig, animate, frameslen(t), init_funcinit, blitTrue, interval50) plt.close() HTML(ani.to_jshtml())这段代码会生成一个双面板动画左面板展示系统状态在相平面中的移动右面板展示Lyapunov函数值随时间的变化通过动画我们可以清晰地看到系统状态如何沿着Lyapunov函数的梯度下降V(x)值如何随时间单调递减最终系统如何收敛到平衡点(V(x)0)6. 实际应用案例倒立摆控制让我们看一个更实际的例子——倒立摆的稳定性分析。倒立摆是一个经典的控制问题我们可以尝试为其设计Lyapunov函数。首先定义倒立摆的动力学方程def inverted_pendulum(x, t, m1.0, l1.0, b0.1, g9.81): theta, theta_dot x dtheta_dt theta_dot dtheta_dot_dt (m*g*l*np.sin(theta) - b*theta_dot) / (m*l**2) return [dtheta_dt, dtheta_dot_dt] # 设计Lyapunov函数候选 def V_pendulum(theta, theta_dot): return 0.5*theta_dot**2 g/l*(1 - np.cos(theta)) # 计算轨迹 t_pendulum np.linspace(0, 10, 1000) initial_angles [np.pi/6, np.pi/4, np.pi/3] # 小角度初始条件 trajectories_pendulum [odeint(inverted_pendulum, [angle, 0], t_pendulum) for angle in initial_angles] # 可视化 plt.figure(figsize(14, 6)) # 相平面 plt.subplot(1, 2, 1) theta np.linspace(-np.pi/2, np.pi/2, 50) theta_dot np.linspace(-2, 2, 50) Theta, Theta_dot np.meshgrid(theta, theta_dot) # 计算Lyapunov函数值 V_values V_pendulum(Theta, Theta_dot) contour plt.contour(Theta, Theta_dot, V_values, levels15, cmapviridis) plt.colorbar(contour, labelV(x)) # 绘制轨迹 colors [r, b, g] for traj, color in zip(trajectories_pendulum, colors): plt.plot(traj[:, 0], traj[:, 1], color, linewidth2) plt.scatter(traj[0, 0], traj[0, 1], colorcolor, s50) plt.xlabel(角度(theta)) plt.ylabel(角速度(theta_dot)) plt.title(倒立摆相平面与Lyapunov函数) plt.grid(True) # 能量变化 plt.subplot(1, 2, 2) for traj, color in zip(trajectories_pendulum, colors): V_traj V_pendulum(traj[:, 0], traj[:, 1]) plt.plot(t_pendulum, V_traj, color, linewidth2, labelf初始角度{np.rad2deg(traj[0, 0]):.1f}°) plt.xlabel(时间(s)) plt.ylabel(V(x)) plt.title(Lyapunov函数随时间变化) plt.legend() plt.grid(True) plt.tight_layout() plt.show()这个例子展示了如何为机械系统设计Lyapunov函数这里使用了总能量在小角度范围内系统是稳定的V(x)递减初始角度越大Lyapunov函数变化越明显提示在实际控制系统中我们通常会设计控制器来确保dV/dt 0。这可以通过在Lyapunov函数分析中引入控制输入u来实现。