From cf58a3b4be80dc86932bd23357dd84ae64cddd19 Mon Sep 17 00:00:00 2001 From: Russ Tedrake Date: Wed, 9 Sep 2026 07:51:09 -0400 Subject: [PATCH] Avoid symbolic expansion when building approximate DP SOS programs --- book/lyapunov/approximate_dp.ipynb | 64 +++++++++++++----------------- 1 file changed, 28 insertions(+), 36 deletions(-) diff --git a/book/lyapunov/approximate_dp.ipynb b/book/lyapunov/approximate_dp.ipynb index 8d82c400..77f25681 100644 --- a/book/lyapunov/approximate_dp.ipynb +++ b/book/lyapunov/approximate_dp.ipynb @@ -311,36 +311,32 @@ " J_dot = J_expr.Jacobian(z).dot(f(z, u))\n", " LHS = J_dot + l(z, u)\n", "\n", + " # Keep the multipliers as Polynomials to avoid expanding their symbolic\n", + " # coefficients when constructing the SOS constraints.\n", " # S procedure for s^2 + c^2 = 1.\n", - " lam_r = prog.NewFreePolynomial(Variables(z), deg).ToExpression()\n", - " S_r = lam_r * (z[0] ** 2 + z[1] ** 2 - 1)\n", + " lam_r = prog.NewFreePolynomial(Variables(z), deg)\n", + " S_r = lam_r * Polynomial(z[0] ** 2 + z[1] ** 2 - 1)\n", " S_Jdot = 0\n", " for i in range(nz):\n", - " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[\n", - " 0\n", - " ].ToExpression()\n", - " S_Jdot += lam * (z[i] - z_max[i]) * (z[i] - z_min[i])\n", + " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[0]\n", + " S_Jdot += lam * Polynomial((z[i] - z_max[i]) * (z[i] - z_min[i]))\n", "\n", " # Enforce Input constraint\n", " u_min = -u_max\n", " for i in range(nu):\n", - " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[\n", - " 0\n", - " ].ToExpression()\n", - " S_Jdot += lam * (u[i] - u_max[i]) * (u[i] - u_min[i])\n", + " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[0]\n", + " S_Jdot += lam * Polynomial((u[i] - u_max[i]) * (u[i] - u_min[i]))\n", " # Enforce Bellman inequality.\n", - " prog.AddSosConstraint(LHS + S_r + S_Jdot)\n", + " prog.AddSosConstraint(Polynomial(LHS, prog.indeterminates()) + S_r + S_Jdot)\n", "\n", - " lam_r = prog.NewFreePolynomial(Variables(z), deg).ToExpression()\n", - " S_r = lam_r * (z[0] ** 2 + z[1] ** 2 - 1)\n", + " lam_r = prog.NewFreePolynomial(Variables(z), deg)\n", + " S_r = lam_r * Polynomial(z[0] ** 2 + z[1] ** 2 - 1)\n", " S_J = 0\n", " for i in range(nz):\n", - " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[\n", - " 0\n", - " ].ToExpression()\n", - " S_J += lam * (z[i] - z_max[i]) * (z[i] - z_min[i])\n", + " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[0]\n", + " S_J += lam * Polynomial((z[i] - z_max[i]) * (z[i] - z_min[i]))\n", " # Enforce that value function is PD\n", - " prog.AddSosConstraint(J_expr + S_r + S_J)\n", + " prog.AddSosConstraint(J + S_r + S_J)\n", "\n", " # J(z0) = 0.\n", " J0 = J_expr.EvaluatePartial(dict(zip(z, z0)))\n", @@ -619,37 +615,33 @@ " LHS = J_dot + l_cost(z, u) * denominator\n", "\n", " lam_deg = Polynomial(LHS).TotalDegree()\n", + " # Keep the multipliers as Polynomials to avoid expanding their symbolic\n", + " # coefficients when constructing the SOS constraints.\n", " # S procedure for s^2 + c^2 = 1.\n", - " lam = prog.NewFreePolynomial(Variables(z), lam_deg).ToExpression()\n", - " S_procedure = lam * (z[1] ** 2 + z[2] ** 2 - 1)\n", + " lam = prog.NewFreePolynomial(Variables(z), lam_deg)\n", + " S_procedure = lam * Polynomial(z[1] ** 2 + z[2] ** 2 - 1)\n", " S_Jdot = 0\n", " for i in np.arange(nz):\n", - " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(lam_deg / 2) * 2))[\n", - " 0\n", - " ].ToExpression()\n", - " S_Jdot += lam * (z[i] - z_max[i]) * (z[i] - z_min[i])\n", + " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(lam_deg / 2) * 2))[0]\n", + " S_Jdot += lam * Polynomial((z[i] - z_max[i]) * (z[i] - z_min[i]))\n", "\n", " # Enforce Input constraint\n", " u_min = -u_max\n", " for i in range(nu):\n", - " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[\n", - " 0\n", - " ].ToExpression()\n", - " S_Jdot += lam * (u[i] - u_max[i]) * (u[i] - u_min[i])\n", + " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[0]\n", + " S_Jdot += lam * Polynomial((u[i] - u_max[i]) * (u[i] - u_min[i]))\n", "\n", - " prog.AddSosConstraint(LHS + S_procedure + S_Jdot)\n", + " prog.AddSosConstraint(Polynomial(LHS, prog.indeterminates()) + S_procedure + S_Jdot)\n", "\n", " # Enforce that value function is PD\n", " S_J = 0\n", - " lam_r = prog.NewFreePolynomial(Variables(z), deg).ToExpression()\n", - " S_r = lam_r * (z[1] ** 2 + z[2] ** 2 - 1)\n", + " lam_r = prog.NewFreePolynomial(Variables(z), deg)\n", + " S_r = lam_r * Polynomial(z[1] ** 2 + z[2] ** 2 - 1)\n", " for i in np.arange(nz):\n", - " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[\n", - " 0\n", - " ].ToExpression()\n", - " S_J += lam * (z[i] - z_max[i]) * (z[i] - z_min[i])\n", + " lam = prog.NewSosPolynomial(Variables(z), int(np.ceil(deg / 2) * 2))[0]\n", + " S_J += lam * Polynomial((z[i] - z_max[i]) * (z[i] - z_min[i]))\n", " # Enforce that value function is PD\n", - " prog.AddSosConstraint(J_expr + S_J + S_r)\n", + " prog.AddSosConstraint(J + S_J + S_r)\n", "\n", " # J(z0) = 0.\n", " J0 = J_expr.EvaluatePartial(dict(zip(z, z0)))\n",