Open jgillis opened 11 months ago
With #3183 implemented, this looks a lot better already. Let's keep this issue open and steer it new something new:
from casadi import *
N = 10
# Physical constants
g = 9.81 # gravitation [m/s^2]
L = 0.2 # pendulum length [m]
m = 1 # pendulum mass [kg]
mcart = 0.5 # cart mass [kg]
T = 2.0 # control horizon [s]
dt = T/N # length of 1 control interval [s]
# System is composed of 4 states
nx = 4
# Construct a CasADi function for the ODE right-hand side
x = MX.sym('x',nx) # states: pos [m], theta [rad], dpos [m/s], dtheta [rad/s]
u = MX.sym('u') # control force [N]
ddpos = ((u+m*L*x[3]*x[3]*sin(x[1])-m*g*sin(x[1])*cos(x[1]))/(mcart+m-m*cos(x[1])*cos(x[1])))
rhs = vertcat(x[2],x[3],ddpos,g/L*sin(x[1])-cos(x[1])*ddpos)
# Continuous system dynamics as a CasADi Function
f = Function('f', [x, u],[rhs])
f = f.expand()
k1 = f(x, u)
k2 = f(x + dt/2 * f(x, u), u)
k3 = f(x + dt/2 * f(x + dt/2 * f(x, u), u), u)
k4 = f(x + dt * f(x + dt/2 * f(x + dt/2 * f(x, u), u), u), u)
x_next = x+dt/6*(k1 +2*k2 +2*k3 +k4)
y = cse(x_next)
g = Function('g',[x,u],[y])
f_recovered =g.find_functions()[0]
f_recovered.disp(True) # Missed cse opportunities