Document author:  itabrah2@illinois.edu

Prompt: I want to create a simple interactive notebook to show students the relation between Q, R and N matrices and pendulum dynamics where the objective is to stabilize an inverted pendulum using LQR control. I want to do this in Julia using Pluto notebooks and the Pluto UI sliders. Please help. The notebook should have three sections -- (a) derivations of EOM using the Symbolics.jl package followed by linearization, (b) implementation of LQR control and (c) interactive part to vary LQR gains and see change in system response.

System: Google Gemini 2.5 Pro

Response

Of course. Here is the complete code for a Pluto notebook that interactively demonstrates LQR control for an inverted pendulum.

You can copy and paste the entire code block below into a single cell in a new Julia Pluto notebook. Pluto's reactive nature will automatically separate it into multiple cells and create the interactive demonstration.^[a]

Julia

### A Pluto.jl notebook for a lecture on inverted pendulum dynamics and LQR control
#
# Sections:
# (a) Derivations of EOM using the Symbolics.jl package followed by linearization
# (b) Implementation of LQR control
# (c) Interactive part to vary LQR gains and see change in system response
#
# This notebook was generated by Gemini.
###

# ------------------------------------------------------------------------------------------
#
# ## 📦 Package Installation and Loading
#
# The first step is to load all the required packages. If they are not installed,
# Julia's package manager will handle the installation automatically.
#
# ------------------------------------------------------------------------------------------

begin
        # Import required packages
        using PlutoUI
        using Symbolics
        using ControlSystems
        using DifferentialEquations
        using LinearAlgebra
        using Plots

        # Set plotting backend
        gr()
        # Set a nice theme for plots
        theme(:dark)
end

# ------------------------------------------------------------------------------------------
#
# # (a) 📜 Derivation of Equations of Motion
#
# In this section, we derive the system's equations of motion (EOM) using Lagrangian mechanics. The [`Symbolics.jl`](https://symbolics.juliasymbols.org/stable/) package allows us to perform these derivations automatically.
#
# The system consists of a cart of mass $M$ and a pendulum of mass $m$ with length $L$. The states are the cart's position $x$ and the pendulum's angle $\theta$ (where $\theta=0$ is the upright position). A force $u$ is applied to the cart.
#
# ------------------------------------------------------------------------------------------

md"""
### Symbolic Variables

First, we define all symbolic variables for our system parameters, states, and their time derivatives.
"""

begin
        @variables t                           # time
        @variables m M L g                     # system parameters
        @variables x(t) θ(t) u(t)              # states and input
        
        # Define differential operator
        D = Differential(t)
        
        # Display variables nicely in LaTeX
        (t, m, M, L, g, x, θ, u, D)
end

md"""
### Lagrangian Mechanics

We use the principle of Lagrangian mechanics to find the EOM. The Lagrangian $\mathcal{L}$ is defined as the difference between the system's kinetic energy ($T$) and potential energy ($V$).

$\mathcal{L} = T - V$

1.  **Kinetic Energy ($T$)**: This is the sum of the kinetic energy of the cart and the pendulum bob.
   * $T_{cart} = \frac{1}{2} M (\dot{x})^2$
   * $T_{pendulum} = \frac{1}{2} m (v_{px}^2 + v_{py}^2)$
   where $x_p = x + L \sin(\theta)$ and $y_p = L \cos(\theta)$.^[b]

2.  **Potential Energy ($V$)**: This is determined by the height of the pendulum bob.
   * $V = m g L \cos(\theta)$

The equations of motion are then found using the Euler-Lagrange equation for each generalized coordinate ($q_i \in \{x, \theta\}$), including the external force $u$ applied to $x$:

$\frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{q_i}} - \frac{\partial \mathcal{L}}{\partial q_i} = F_{q_i}$
"""

begin
        # Position of pendulum bob
        x_p = x + L * sin(θ)
        y_p = L * cos(θ)

        # Velocity of pendulum bob
        v_px = D(x_p)
        v_py = D(y_p)

        # 1. Kinetic Energy (T)
        T_cart = 0.5 * M * D(x)^2
        T_pend = 0.5 * m * (v_px^2 + v_py^2)
        T = T_cart + T_pend

        # 2. Potential Energy (V)
        V = m * g * L * cos(θ)

        # 3. Lagrangian (L)
        Lag = T - V

        # 4. Euler-Lagrange Equations
        # Generalized coordinates q = [x, θ]
        # Forces F = [u, 0]
        EL_x = D(gradient(Lag, [D(x)])[1]) - gradient(Lag, [x])[1] ~ u
        EL_θ = D(gradient(Lag, [D(θ)])[1]) - gradient(Lag, [θ])[1] ~ 0

        # Display the two EOMs
        (EL_x, EL_θ)
end

md"""
### Solving for Accelerations

The equations above are implicit. To get our standard state-space form, we need to solve for the accelerations $\ddot{x}$ and $\ddot{\theta}$.
"""

begin
        # Solve for the second derivatives (accelerations)
        sol = Symbolics.solve_for^[c]([EL_x, EL_θ^[d]], [D(D(x)), D(D(θ))])
        
        x_dd_expr = sol[1]
        θ_dd_expr = sol[2]

        md"""
        The solved expression for cart acceleration $\ddot{x}$ is:
        $$ \ddot{x} = $(latexify(x_dd_expr)) $$

        The solved expression for pendulum angular acceleration $\ddot{\theta}$ is:
        $$ \ddot{\theta} = $(latexify(θ_dd_expr)) $$
        """
end

md"""
### System Linearization

The LQR controller requires a linear time-invariant (LTI) model of the form $\dot{\mathbf{z}} = A\mathbf{z} + B\mathbf{u}$. We linearize our nonlinear system around the unstable equilibrium point: pendulum perfectly upright and stationary at the origin.

The state vector is $\mathbf{z} = [x, \dot{x}, \theta, \dot{\theta}]^T$.
The equilibrium point is $\mathbf{z}_{eq} = [0, 0, 0, 0]^T$ and $u_{eq} = 0$.

The matrices $A$ and $B$ are the Jacobians of the system dynamics with respect to the state and input, evaluated at the equilibrium point.

$A = \frac{\partial f}{\partial \mathbf{z}} \bigg|_{\mathbf{z}_{eq}, u_{eq}}$ and $B = \frac{\partial f}{\partial u} \bigg|_{\mathbf{z}_{eq}, u_{eq}}$
"""

begin
        # Define the symbolic state vector and nonlinear dynamics function f
        z = [x, D(x), θ, D(θ)]
        f_sym = [
                D(x),
                x_dd_expr,
                D(θ),
                θ_dd_expr
        ]
        
        # Compute the Jacobians symbolically
        A_sym = Symbolics.jacobian(f_sym, z)
        B_sym = Symbolics.jacobian(f_sym, [u])

        # Define the equilibrium point
        eq_point = Dict(x => 0, D(x) => 0, θ => 0, D(θ) => 0, u => 0)

        # Substitute the equilibrium point values into the Jacobians
        A_eq = substitute(A_sym, eq_point)
        B_eq = substitute(B_sym, eq_point)

        md"""
        The symbolic linearized state matrix $A$ is:
        $$ A = $(latexify(A_eq)) $$

        The symbolic linearized input matrix $B$ is:
        $$ B = $(latexify(B_eq)) $$
        """
end

# ------------------------------------------------------------------------------------------
#
# # (b) 🤖 LQR Controller Implementation
#
# The Linear-Quadratic Regulator (LQR) is an optimal control strategy. It finds the feedback control law $u = -K\mathbf{z}$ that minimizes a quadratic cost function:
#
# $J = \int_{0}^{\infty} (\mathbf{z}^T Q \mathbf{z} + u^T R u) \,dt$
#
# * **Q matrix**: A positive semi-definite matrix that penalizes deviations in the state variables from the equilibrium.
# * **R matrix**: A positive definite matrix that penalizes control effort.
#
# The trade-off is key: large values in $Q$ lead to aggressive control, while large values in $R$ lead to conservative, energy-saving control.
#
# ------------------------------------------------------------------------------------------

md"""
### Numerical System Model

Let's define numerical values for our system parameters and compute the numerical $A$ and $B$ matrices.
"""

begin
        # Define a dictionary with numerical parameters
        # We use `Float64` to ensure the resulting matrices are not symbolic
        params = Dict(
                M => 1.0,      # Cart mass (kg)
                m => 0.2,      # Pendulum mass (kg)
                L => 0.5,      # Pendulum length (m)
                g => 9.81      # Gravitational acceleration (m/s²)
        )

        # Substitute numerical values into the symbolic A and B matrices
        A_num = Float64.(substitute(A_eq^[e], params))
        B_num = Float64.(substitute(B_eq^[f], params))
        
        # Create the state-space model using ControlSystems.jl
        sys_linear = ss(A_num, B_num, I(4), 0)
end

md"""
Using the parameters $M = $(params[M])\,$kg, $m = $(params[m])\,$kg, $L = $(params[L])\,$m, and $g = $(params[g])\,$m/s², our linearized system matrices are:

$A = $(sprint(show, "text/plain", A_num))$

$B = $(sprint(show, "text/plain", B_num))$

The system is unstable, which we can confirm by checking the eigenvalues of the $A$ matrix. A positive eigenvalue indicates instability.

Eigenvalues of A: **$(round.(eigvals(A_num), digits=2))**
"""


# ------------------------------------------------------------------------------------------
#
# # (c) 🕹️ Interactive LQR Simulation
#
# Now for the fun part! Use the sliders below to adjust the LQR cost matrices $Q$ and $R$. Observe how the system's response changes.
#
# * **Position Cost ($Q_{11}$)**: How much to penalize the cart for not being at the origin ($x=0$).
# * **Angle Cost ($Q_{33}$)**: How much to penalize the pendulum for not being upright ($\theta=0$).
# * **Control Cost ($R_{11}$)**: How much to penalize using the motor force ($u$).
#
# We will simulate the response of both the **linearized model** and the **original nonlinear model** to the same LQR controller to see how they compare.
#
# ------------------------------------------------------------------------------------------

md"""
### LQR Cost Sliders
"""

@bind q1_slider Slider(0.1:0.1:20, default=1.0, show_value=true)
@bind q3_slider Slider(0.1:0.1:20, default=10.0, show_value=true)
@bind r_slider Slider(0.01:0.01:2.0, default=0.1, show_value=true)

md"""
**Position Cost ($Q_{11}$)**: $( PlutoUI.Show(q1_slider) )
**Angle Cost ($Q_{33}$)**: $( PlutoUI.Show(q3_slider) )
**Control Cost ($R$)**: $( PlutoUI.Show(r_slider) )
"""

begin
        # Define Q and R matrices based on slider values
        # We penalize position (z1) and angle (z3). Velocity penalties (z2, z4) are often less critical for stabilization.
        Q = diagm([q1_slider, 0, q3_slider, 0])
        R = [r_slider]

        # Compute the LQR gain K for the linear system
        # The lqr function solves the algebraic Riccati equation to find K
        K, _, _ = lqr(sys_linear, Q, R)
        
        # Display the calculated gain matrix
        md"""
        The resulting LQR gain matrix $K$ is: **$(round.(K, digits=2))**
        """
end

md"""
### System Simulation

We now simulate the closed-loop system dynamics for an initial condition where the pendulum is slightly tilted.

* Initial state: $\mathbf{z}_0 = [x=0, \dot{x}=0, \theta=0.1, \dot{\theta}=0]^T$.
* Control law: $u = -K \mathbf{z}$
"""

begin
        # ---- 1. Simulation Setup ----
        x₀ = 0.0    # initial cart position
        θ₀ = 0.1    # initial pendulum angle (in radians, approx 5.7 degrees)
        z₀ = [x₀, 0.0, θ₀, 0.0]  # Initial state vector
        t_span = (0.0, 10.0)    # Simulation time span

        # ---- 2. Nonlinear System Dynamics ----
        # Here we build the function for the ODE solver using the symbolic expressions
        # This avoids manual transcription errors and is highly reusable.
        # The `build_function` call generates highly optimized code.
        f_nonlinear_oop = build_function(
                f_sym,
                [z...], [u], [M, m, L, g],
                t,
                expression=Val{false}
        )[1]

        # Define the closed-loop function for the nonlinear system
        function nonlinear_dynamics!(dz, z, p, t)
                # Unpack parameters and gain matrix
                M_val, m_val, L_val, g_val, K_mat = p
                
                # Calculate control input
                u_val = -K_mat * z
                
                # Call the generated function
                # Note: The generated function expects a scalar u, so we extract it.
                dz .= f_nonlinear_oop(z, u_val[1], [M_val, m_val, L_val, g_val], t)
        end
        
        # ---- 3. Linear System Dynamics ----
        # Define the closed-loop function for the linearized system
        function linear_dynamics!(dz, z, p, t)
                A_mat, B_mat, K_mat = p
                dz .= (A_mat - B_mat * K_mat) * z
        end
        
        # ---- 4. Solve the ODEs ----
        # Solve nonlinear system
        p_nonlin = (params[M], params[m], params[L], params[g], K)
        prob_nonlin = ODEProblem(nonlinear_dynamics!, z₀, t_span, p_nonlin)
        sol_nonlin = solve(prob_nonlin, Tsit5(), saveat=0.05)

        # Solve linear system
        p_lin = (A_num, B_num, K)
        prob_lin = ODEProblem(linear_dynamics!, z₀, t_span, p_lin)
        sol_lin = solve(prob_lin, Tsit5(), saveat=0.05)
end

md"""
### 📊 Simulation Results

The plots below show the cart position, pendulum angle, and control force over time. The solid line represents the **nonlinear model**, while the dashed line represents the **linear model**.
"""

begin
        # Plot 1: Cart Position (x)
        p1 = plot(
^[g]                sol_lin.t, [s[1] for s in sol_lin.u],
                label="Linear", linestyle=:dash, lw=2,
                title="Cart Position", xlabel="Time (s)", ylabel="Position (m)"
        )
        plot!(
                p1, sol_nonlin.t, [s[1] for s in sol_nonlin.u],
                label="Nonlinear", lw=2
        )

        # Plot 2: Pendulum Angle (θ)
        p2 = plot(
                sol_lin.t, [rad2deg(s[3]) for s in sol_lin.u],
                label="Linear", linestyle=:dash, lw=2,
                title="Pendulum Angle", xlabel="Time (s)", ylabel="Angle (°)"
        )
        plot!(
                p2, sol_nonlin.t, [rad2deg(s[3]) for s in sol_nonlin.u],
                label="Nonlinear", lw=2
        )

        # Plot 3: Control Input (u)
        # u(t) = -K * z(t)
        u_lin = [-K * s for s in sol_lin.u]
        u_nonlin = [-K * s for s in sol_nonlin.u]
        p3 = plot(
                sol_lin.t, u_lin,
                label="Linear", linestyle=:dash, lw=2,
                title="Control Input", xlabel="Time (s)", ylabel="Force (N)"
        )
        plot!(
                p3, sol_nonlin.t, u_nonlin,
                label="Nonlinear", lw=2
        )
        
        # Combine plots into one layout
        plot(p1, p2, p3, layout=(3,1), size=(800, 700), legend=:topright)
end

md"""
### 🤔 Interpretation

* **High Angle Cost ($Q_{33}$)**: If you increase the penalty on the angle, the controller will act more aggressively to bring $\theta$ to zero, often resulting in larger control inputs and more overshoot in the cart's position.
* **High Position Cost ($Q_{11}$)**: Increasing the penalty on position forces the cart to return to the origin more quickly.
* **High Control Cost ($R$)**: If you increase the penalty on the control force, the controller becomes "lazier" to save energy. This results in a slower response and smaller control signals. If $R$ is too high, the controller may not be strong enough to stabilize the pendulum at all.
* **Linear vs. Nonlinear**: For small initial angles, the linear and nonlinear responses are very similar. If you were to increase the initial angle `θ₀` in the code, you would see the linear model become a less accurate predictor of the true nonlinear behavior.
"""