使用 drake 进行吸引力区域分析

using drake for region of attraction analysis

我正在使用 pydrake 进行吸引区域 (RoA) 分析。有一种涉及三角函数的非线性系统让我陷入困境。 post编辑了一个类似的问题

正如 Russ 在 post 中指出的那样,有 2 个解决方案:1) 转移到另一个状态-space; 2) 使用泰勒展开。我正在将第二种方法应用于钟摆箱以获得感觉。系统动力学在 here 中给出。这是我拥有的:

  1. 将系统线性化,找到一个LQR,其cost-to-go函数将有助于提出李亚普诺夫候选:
# Prepare
builder = DiagramBuilder()
pendulum = builder.AddSystem(PendulumPlant())
context = pendulum.CreateDefaultContext()
pendulum.get_input_port(0).FixValue(context, [0])
context.SetContinuousState([np.pi, 0])

# linearize the system at [pi, 0]
linearized_pendulum = Linearize(pendulum, context)

# LQR
Q = np.diag((10., 1.))
R = [1.]
(K, S) = LinearQuadraticRegulator(linearized_pendulum.A(), linearized_pendulum.B(), Q, R)
  1. 使用 LQR 设计的控制输入设置摆的非线性动力学
# constants
pi = np.pi
g = 9.81
m = 1.
l = 0.5
d = 0.1

# non-linear dyamics
prog = MathematicalProgram()
x = prog.NewIndeterminates(2, "x")
x1 = x[0]
x2 = x[1]
Tsin = -(x1-pi) + (x1-pi)**3/6  # 3rd order taylor expansion of sin(x) around the fix point
r = [pi, 0]  # reference
u = -K.dot(x)[0] # control input
u2 = -K.dot(x-r)[0] # control input with reference
fn = [x2, (u-d*x2-Tsin*m*g*l)/(m*l*l)]
  1. 用s-procedure进行RoA分析,修正一个Lyapunov函数,求拉格朗日乘数
# cost-to-go of LQR as Lyapunov candidate
V = x.dot(S.dot(x)) 
Vdot = Jacobian([V], x).dot(fn)[0]

# Define the Lagrange multiplier.
lambda_ = prog.NewSosPolynomial(Variables(x), 2)[0].ToExpression()

# Optimization setup
rho = 2.
prog.AddSosConstraint(-Vdot + lambda_*(V-rho))
prog.AddSosConstraint(V)
prog.AddSosConstraint(lambda_)

# Print result
result = Solve(prog)
if (result.is_success()):
    print(f"Verified that {str(V)} < {rho} is in the region of attraction.")
    #print(Polynomial(result.GetSolution(lambda_)))
else:
    print("failed")

无论我为 $rho$ 设置的有多小,s-procedure 总是 returns 失败。但是,从here可知,这个例子的RoA至少为2。

我不知道这是一个数值问题还是我的表述有误。如果您能对此有所了解,将不胜感激。谢谢!

您的候选李亚普诺夫函数不正确。目前您使用

V = x.dot(S.dot(x)) 

但我认为你应该使用

V = (x - r).dot(S.dot(x - r))

顺便说一句,作为一个轶事(所以对此持保留态度)。当我使用方法 1 处理一个更复杂的问题时,我取得了成功(即创建一个新变量 xbar = x - r,然后使用 xbar 重写整个动力学,并将 xbar 视为不确定项)。在我尝试方法 2 的问题中,它给了我数字错误。因此,如果方法 2 对您的问题不起作用,您可以尝试方法 1。