Stabilize an Inverted Pendulum with LQR

Simulate a cart-pole system and design an LQR controller to balance it.

Introduction

30 Second Summary

Think about balancing a broomstick on your palm. It falls the instant you stop paying attention. Every self-balancing robot, rocket landing system, and Segway faces the same challenge.

In this project, you will build a complete control pipeline in a Jupyter notebook. You will simulate an inverted pendulum, watch it fall, prove mathematically why it falls, and then design a controller that catches it.

What You'll Build

You will run a simulation, watch a pendulum crash within two seconds, then activate your LQR controller and watch the same pendulum balance itself upright indefinitely.

By the end of this project, you'll have:

  • A nonlinear simulation you can run to watch the uncontrolled pendulum diverge past 90 degrees and crash, plotted as angle vs time with cart position drifting off-screen.
  • A phase portrait and eigenvalue analysis that prove mathematically why the system is unstable, showing trajectories spiraling away from the upright equilibrium.
  • A working LQR feedback controller that stabilizes the pendulum at its unstable upright equilibrium, with overlay plots comparing the controlled and uncontrolled responses side by side.
  • Secret Mission: Disturbance rejection testing where you throw random impulse kicks at the balanced pendulum and tune the Q/R cost matrices to explore the trade-off between aggressive and smooth control.

Are there any prerequisites?

This project assumes basic Python proficiency, comfort with linear algebra (matrices and eigenvalues), and familiarity with differential equations. You do not need prior experience with control theory.

Before We Start

Before we dive in, let's lock in what you're about to build and why it matters for your robotics learning.

Set Up Your Environment

Control theory is a computational discipline. You will write physics equations, solve differential equations numerically, and plot the results to build intuition. That workflow demands a proper scientific computing environment.

In this step, you will create an isolated Python environment and install the same stack used by Steve Brunton in his Control Bootcamp: NumPy for linear algebra, SciPy for differential equation solvers and the Riccati solver, and Matplotlib for visualization.

In this step, get ready to:
  • Create and activate a virtual environment.
  • Install NumPy, SciPy, Matplotlib, and Jupyter Notebook.
  • Launch Jupyter and verify all imports work.
Create a virtual environment

A virtual environment keeps your project's packages isolated from the rest of your system. This means your control theory libraries will not interfere with other Python projects on your machine.

  • Open your terminal. Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Terminal and press Enter to open it.
  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your virtual environment and activate it by running these commands:
python3 -m venv control-lab
source control-lab/bin/activate

Your terminal prompt now shows (control-lab) at the beginning. This confirms you are working inside the isolated environment.

What do these commands do?

The first command creates a folder called control-lab on your Desktop containing a fresh Python installation with its own package directory.

The second command activates it so that any pip install commands install packages into this environment only.

Don't see (control-lab) in your prompt?

Make sure you ran both commands in the same terminal window. If you closed the terminal between them, the activation is lost.

Try running source ~/Desktop/control-lab/bin/activate using the full path.

Still stuck? Why isn't my virtual environment activating?

Install the scientific Python stack

With your environment active, you can now install the four packages this project needs. Jupyter Notebook gives you the interactive cell-by-cell workflow that makes iterating on simulations fast.

  • Install all four packages by running this command:
pip install numpy scipy matplotlib notebook

This will take a minute or two. You will see a stream of download and installation progress messages.

  • Wait for the install to finish. Confirm you see a Successfully installed message listing numpy, scipy, matplotlib, and notebook.

Seeing a permissions error or build failure?

If you see a permissions error, check that (control-lab) appears in your prompt. Installing without an active virtual environment requires admin permissions.

If SciPy fails to build, confirm your Python version is 3.12 or higher by running python3 --version. SciPy 1.18.0 requires Python 3.12+.

Need more help? My pip install is failing for numpy scipy matplotlib notebook.

Launch Jupyter and verify your imports

Everything is installed. Now you will launch Jupyter Notebook, create your project notebook, and run a test cell to confirm the entire stack works together.

  • Launch Jupyter Notebook by running this command:
jupyter notebook

What happens when you run this?

This command starts a local web server and opens your default browser to the Jupyter file browser. Your terminal will be occupied running the server, so leave it open.

The file browser shows your Desktop directory. You will create your notebook here.

  • In the Jupyter file browser in your browser, click New in the top-right corner.
  • Select Python 3 to create a new notebook.
  • Click the title at the top of the notebook (it says Untitled) and rename it to cart_pole_control.
  • In the first empty cell, paste the following code:
# Cell 1: Imports and Setup
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_are
import matplotlib.pyplot as plt
print("Setup complete!")

What does this code do?

  • numpy as np gives you fast array math and linear algebra operations.
  • solve_ivp is SciPy's differential equation solver. You will use it to simulate the cart-pole physics.
  • solve_continuous_are solves the Riccati equation at the heart of LQR control design.
  • matplotlib.pyplot as plt handles all the plotting: time series, phase portraits, and comparison figures.

Before you run this cell, what do you expect to see if all four packages installed correctly?

  • Run the cell by pressing Shift+Enter.

You should see Setup complete! printed below the cell with no error messages. This confirms NumPy, SciPy, and Matplotlib are all importable and ready to use.

Seeing a ModuleNotFoundError?

This means the notebook is using a different Python than the one in your virtual environment. Close the Jupyter browser tab and stop the server in your terminal with Ctrl+C.

Confirm (control-lab) is still in your prompt, then run jupyter notebook again. The server must be launched from inside the activated environment.

Still seeing the error? Help me fix ModuleNotFoundError in my Jupyter notebook.

✔️ Awesome, I've got everything!

Great. Double check you see Setup complete! in the output, then you are ready to move on.

ⓧ I'd like to double check the full code

Your first cell in cart_pole_control.ipynb should contain exactly this:

# Cell 1: Imports and Setup
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_are
import matplotlib.pyplot as plt
print("Setup complete!")

Your scientific computing environment is locked and loaded. Next up, you will write the full nonlinear physics of the cart-pole system and watch what happens when you let the pendulum go.

Simulate the Uncontrolled Cart-Pole

Before you can control a system, you need to see it fail on its own. A pendulum balanced upright is inherently unstable. Even a tiny nudge sends it crashing down.

In this step, you will write the full nonlinear equations of motion for the cart-pole system and simulate what happens when the pendulum starts slightly off-vertical. This gives you the "problem" that motivates everything that follows.

In this step, get ready to:
  • Define the physical parameters and nonlinear equations of motion for the cart-pole system.
  • Simulate 3 seconds of uncontrolled motion starting from a tiny tilt.
  • Plot the pendulum angle and cart position to see the system fail.
Define the system parameters

The cart-pole system has four physical constants: the mass of the cart, the mass of the pendulum bob, the length of the pendulum rod, and gravity. These values determine how the system behaves.

  • In your notebook, create a new cell below Cell 1. Define the system parameters by pasting this code:
# Cell 2: System Parameters
M = 1.0      # cart mass (kg)
m = 0.1      # pendulum mass (kg)
L = 0.5      # pendulum length (m)
g = 9.81     # gravity (m/s^2)

What do these values represent?

The cart weighs 1 kg and rolls horizontally. The pendulum bob (0.1 kg) sits at the end of a 0.5 m rod attached to the cart. Gravity pulls at 9.81 m/s². The pendulum is ten times lighter than the cart, which makes it easy to disturb but hard to catch.

  • Run the cell by pressing Shift+Enter. The cell executes silently (no output) and your cursor moves to the next cell.
Write the nonlinear equations of motion

The cart-pole dynamics come from the Euler-Lagrange equations. These describe how forces and accelerations couple between the cart and pendulum through the shared pivot point.

  • Create a new cell below Cell 2. Define the dynamics function by pasting this code:
# Cell 3: Nonlinear Equations of Motion
def cart_pole_dynamics(t, state, F=0.0):
    """
    Full nonlinear dynamics of the cart-pole system.
    State: [x, x_dot, theta, theta_dot]
    theta = 0 is the UPRIGHT (unstable) equilibrium.
    F is the horizontal force applied to the cart.
    """
    x, x_dot, theta, theta_dot = state
    sin_t = np.sin(theta)
    cos_t = np.cos(theta)

    denom = M + m * sin_t**2

    x_ddot = (F + m * sin_t * (L * theta_dot**2 - g * cos_t)) / denom
    theta_ddot = ((M + m) * g * sin_t - F * cos_t
                  - m * L * theta_dot**2 * sin_t * cos_t) / (L * denom)

    return [x_dot, x_ddot, theta_dot, theta_ddot]

What does this function do?

The function takes the current state [x, x_dot, theta, theta_dot] and returns the derivatives (velocities and accelerations). theta = 0 means the pendulum is perfectly upright.

  • denom captures the nonlinear coupling between cart and pendulum. It changes with the pendulum angle.
  • x_ddot is the cart's horizontal acceleration, driven by the applied force and the pendulum's swinging.
  • theta_ddot is the angular acceleration of the pendulum. Gravity pulls it down, and the cart's motion tugs it sideways.
  • Run the cell with Shift+Enter.
  • Verify the function works by creating a new cell and running this quick test:
cart_pole_dynamics(0, [0, 0, 0.1, 0], 0.0)

You should see a list of four numbers. The third value (angular acceleration) is positive, meaning gravity is already pulling the pendulum away from vertical.

Seeing a NameError?

  • If you see NameError: name 'M' is not defined, go back and run Cell 2 (the parameters cell) first.
  • If you see NameError: name 'np' is not defined, go back and run Cell 1 (the imports cell) first.

Why can't my dynamics function find the parameters?

  • Delete the test cell once you have confirmed the function works. You will not need it again.
Simulate the uncontrolled pendulum

Now you will set the pendulum just 0.1 radians (about 5.7 degrees) off vertical and let go. No force is applied to the cart. The simulation runs for 3 seconds using solve_ivp from SciPy.

  • Create a new cell below Cell 3. Run the uncontrolled simulation by pasting this code:
# Cell 4: Simulate Uncontrolled System
theta0 = 0.1  # initial tilt (radians) - about 5.7 degrees
state0 = [0.0, 0.0, theta0, 0.0]  # [x, x_dot, theta, theta_dot]
t_span = (0, 3)
t_eval = np.linspace(0, 3, 500)

sol_uncontrolled = solve_ivp(
    cart_pole_dynamics, t_span, state0,
    t_eval=t_eval, args=(0.0,), method='RK45'
)

What does this code do?

  • state0 = [0.0, 0.0, theta0, 0.0] sets the initial condition: cart at the origin, not moving, pendulum tilted 0.1 radians from vertical, no angular velocity.
  • solve_ivp numerically integrates the equations of motion forward in time using the RK45 method (a standard fourth-order Runge-Kutta solver).
  • args=(0.0,) passes zero force to the dynamics function. Nobody is pushing the cart.
  • Run the cell with Shift+Enter. The simulation runs in under a second with no output. The result is stored in sol_uncontrolled.
Plot the uncontrolled response

Time to see what happened. You will create a 2x1 subplot: pendulum angle on top, cart position on the bottom.

  • Create a new cell below Cell 4. Visualize the uncontrolled response by pasting this code:
# Cell 5: Plot Uncontrolled Response
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(sol_uncontrolled.t, np.degrees(sol_uncontrolled.y[2]), 'r', linewidth=2)
axes[0].axhline(y=0, color='k', linestyle='--', alpha=0.3)
axes[0].set_ylabel('Pendulum Angle (degrees)')
axes[0].set_title('Uncontrolled Cart-Pole: The Pendulum Falls')
axes[0].grid(True, alpha=0.3)

axes[1].plot(sol_uncontrolled.t, sol_uncontrolled.y[0], 'b', linewidth=2)
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Cart Position (m)')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('uncontrolled_response.png', dpi=150)
plt.show()

What does this plotting code do?

  • np.degrees(sol_uncontrolled.y[2]) converts the pendulum angle from radians to degrees for readability.
  • The dashed horizontal line at zero marks the upright equilibrium. Any deviation from this line means the pendulum is falling.
  • plt.savefig('uncontrolled_response.png', dpi=150) saves the figure to disk so you can reference it later.

Before you run this cell, what do you think will happen to the pendulum angle over 3 seconds? Will the tiny 5.7-degree tilt stay small, oscillate, or grow?

  • Run the cell with Shift+Enter.

You should see a 2-panel figure. The top plot shows the pendulum angle rocketing past 90 degrees within about 1.5 seconds. The bottom plot shows the cart drifting off-screen as the pendulum drags it. A tiny 5.7-degree nudge was enough to destroy the system.

Why does it fall so fast?

The upright position is an unstable equilibrium. Any deviation, no matter how small, gets amplified by gravity. The system does not just drift. It diverges exponentially. You will prove this mathematically with eigenvalue analysis in a later step.

Not seeing any plot?

  • Make sure you ran Cell 4 (the simulation cell) before Cell 5. The variable sol_uncontrolled must exist before you can plot it.
  • If the plot appears but looks empty or flat, check that Cell 2 (parameters) was run. Without g = 9.81, gravity is not defined and the dynamics produce zeros.

Why isn't my uncontrolled simulation plotting anything?

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with all five cells.

ⓧ I'd like to double check the full code

Here is the full notebook so far. Scroll up and confirm each cell matches in order.

Cell 1: Imports and Setup

# Cell 1: Imports and Setup
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_are
import matplotlib.pyplot as plt
print("Setup complete!")

Cell 2: System Parameters

# Cell 2: System Parameters
M = 1.0      # cart mass (kg)
m = 0.1      # pendulum mass (kg)
L = 0.5      # pendulum length (m)
g = 9.81     # gravity (m/s^2)

Cell 3: Nonlinear Equations of Motion

# Cell 3: Nonlinear Equations of Motion
def cart_pole_dynamics(t, state, F=0.0):
    x, x_dot, theta, theta_dot = state
    sin_t = np.sin(theta)
    cos_t = np.cos(theta)
    denom = M + m * sin_t**2
    x_ddot = (F + m * sin_t * (L * theta_dot**2 - g * cos_t)) / denom
    theta_ddot = ((M + m) * g * sin_t - F * cos_t
                  - m * L * theta_dot**2 * sin_t * cos_t) / (L * denom)
    return [x_dot, x_ddot, theta_dot, theta_ddot]

Cell 4: Simulate Uncontrolled System

# Cell 4: Simulate Uncontrolled System
theta0 = 0.1
state0 = [0.0, 0.0, theta0, 0.0]
t_span = (0, 3)
t_eval = np.linspace(0, 3, 500)
sol_uncontrolled = solve_ivp(
    cart_pole_dynamics, t_span, state0,
    t_eval=t_eval, args=(0.0,), method='RK45'
)

Cell 5: Plot Uncontrolled Response

# Cell 5: Plot Uncontrolled Response
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
axes[0].plot(sol_uncontrolled.t, np.degrees(sol_uncontrolled.y[2]), 'r', linewidth=2)
axes[0].axhline(y=0, color='k', linestyle='--', alpha=0.3)
axes[0].set_ylabel('Pendulum Angle (degrees)')
axes[0].set_title('Uncontrolled Cart-Pole: The Pendulum Falls')
axes[0].grid(True, alpha=0.3)
axes[1].plot(sol_uncontrolled.t, sol_uncontrolled.y[0], 'b', linewidth=2)
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Cart Position (m)')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('uncontrolled_response.png', dpi=150)
plt.show()

You have seen the problem. A tiny tilt destroys the system in under two seconds. In the next step, you will visualize the state space to understand WHY trajectories flee from the upright position.

Visualize the State Space

Your time-series plot from the last step shows the pendulum crashing, but it does not explain why. The angle shoots up, the cart drifts, and all you can say is "it fell."

A phase portrait plots angle against angular velocity on a single plane, revealing the geometric structure of the dynamics. You will see the unstable equilibrium as a saddle point where all nearby trajectories are repelled outward. This is the Brunton/Tedrake approach: build geometric intuition from plots, not just equations.

In this step, get ready to:
  • Simulate the system from a grid of 49 initial conditions near the upright position.
  • Plot all trajectories on a single phase portrait with directional arrows.
  • Add a vector field that makes the saddle-point structure visually obvious.
Simulate and plot trajectories from multiple initial conditions

In the last step, you ran one simulation from a single starting angle. To reveal the geometric structure, you need to see how many trajectories behave at once. A 7x7 grid of starting angles and angular velocities gives you 49 trajectories.

  • Create a new cell below Cell 5 in your notebook and paste the following code:
# Cell 6: Phase Portrait - Multiple Trajectories
fig, ax = plt.subplots(1, 1, figsize=(8, 8))

initial_conditions = []
for theta_ic in np.linspace(-0.3, 0.3, 7):
    for theta_dot_ic in np.linspace(-0.5, 0.5, 7):
        initial_conditions.append([0.0, 0.0, theta_ic, theta_dot_ic])

for ic in initial_conditions:
    sol = solve_ivp(
        cart_pole_dynamics, (0, 2), ic,
        t_eval=np.linspace(0, 2, 200), args=(0.0,), method='RK45'
    )
    ax.plot(sol.y[2], sol.y[3], 'r', alpha=0.4, linewidth=0.8)
    mid = len(sol.t) // 3
    if mid + 1 < len(sol.t):
        ax.annotate('', xy=(sol.y[2][mid+1], sol.y[3][mid+1]),
                    xytext=(sol.y[2][mid], sol.y[3][mid]),
                    arrowprops=dict(arrowstyle='->', color='r', lw=1.5))

What does this code do?

This builds a family of 49 trajectories by sweeping initial angles from -0.3 to 0.3 radians and angular velocities from -0.5 to 0.5 rad/s.

  • Each trajectory is simulated for 2 seconds using solve_ivp with no external force.
  • The ax.plot(sol.y[2], sol.y[3]) line plots angle (theta) on the x-axis and angular velocity (theta_dot) on the y-axis. This is the phase portrait.
  • The ax.annotate call places an arrow at the one-third point of each trajectory, showing which direction time flows.
Add the vector field to reveal the dynamics

Trajectories show individual paths, but a vector field shows the direction and strength of the dynamics at every point in the plane. This makes the saddle-point structure impossible to miss.

  • In the same cell, paste the following code directly below what you just added:
theta_range = np.linspace(-0.4, 0.4, 15)
theta_dot_range = np.linspace(-1.0, 1.0, 15)
THETA, THETA_DOT = np.meshgrid(theta_range, theta_dot_range)
U = THETA_DOT
V = np.zeros_like(THETA)
for i in range(THETA.shape[0]):
    for j in range(THETA.shape[1]):
        state_ij = [0, 0, THETA[i, j], THETA_DOT[i, j]]
        derivs = cart_pole_dynamics(0, state_ij, F=0.0)
        V[i, j] = derivs[3]

ax.quiver(THETA, THETA_DOT, U, V, alpha=0.3, color='gray')
ax.plot(0, 0, 'ko', markersize=10, label='Unstable equilibrium')
ax.set_xlabel(r'$\theta$ (rad)')
ax.set_ylabel(r'$\dot{\theta}$ (rad/s)')
ax.set_title('Phase Portrait: Trajectories Flee the Unstable Equilibrium')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('phase_portrait.png', dpi=150)
plt.show()

What does this code do?

  • A 15x15 grid of points spans a wider region than the trajectories. At each point, the code computes the derivative (angular acceleration) by calling cart_pole_dynamics directly.
  • U = THETA_DOT is the horizontal component of the vector (how theta changes). V[i, j] = derivs[3] is the vertical component (how theta_dot changes). Together they form the quiver arrows.
  • The black dot at the origin marks the unstable upright equilibrium. All arrows point away from it, confirming the saddle-point structure.
Run the cell and observe the phase portrait

Before you run this cell, think about what you expect to see. Will the trajectories converge toward the origin, or will they flee from it?

  • Run the cell by pressing Shift+Enter.

You should see a phase portrait with red trajectories radiating away from the origin in all directions. The gray quiver arrows confirm the flow: everything near the origin is pushed outward. This is the saddle-point structure of the unstable equilibrium.

Why is this a saddle point?

A saddle point has directions that attract and directions that repel. For the inverted pendulum, there is no direction that attracts. Every small perturbation grows. The quiver arrows make this viscerally clear in a way that time-series plots cannot.

This is exactly why feedback control is necessary. No passive mechanism can stabilize this equilibrium.

Not seeing the plot?

  • Make sure both code blocks are in the same cell. If you split them across two cells, the second cell will try to use ax before the figure is displayed.
  • Check that Cells 1-3 have already been executed. The code calls cart_pole_dynamics which was defined in Cell 3.
  • If the plot appears but looks empty (no red lines), verify that the initial_conditions list is not empty by adding print(len(initial_conditions)) before the plotting loop. It should print 49.

Why isn't my phase portrait showing trajectories?

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook and that phase_portrait.png appeared in your working directory.

ⓧ I'd like to double check the full code

Your Cell 6 should contain the two blocks you pasted above in sequence. Scroll up and confirm:

  • The cell starts with # Cell 6: Phase Portrait - Multiple Trajectories and the figure/initial-conditions setup.
  • The trajectory simulation loop with solve_ivp and ax.annotate comes next.
  • The quiver/vector field section starts with theta_range = np.linspace(-0.4, 0.4, 15).
  • The cell ends with plt.savefig('phase_portrait.png', dpi=150) and plt.show().

You now have visual proof that the upright equilibrium repels everything around it. Next, you will prove this mathematically by computing the eigenvalues that quantify exactly how fast the system diverges.

Linearize and Prove Instability

You have seen the pendulum fall in simulation and watched trajectories spiral away from the equilibrium in the phase portrait. Now it is time to prove why mathematically.

Linearization gives you the A and B matrices that describe the system near the upright position. The eigenvalues of A reveal the exact instability rate. The controllability matrix answers the crucial question: can feedback actually fix this?

In this step, get ready to:
  • Derive the linearized A and B matrices at the upright equilibrium.
  • Compute eigenvalues to prove the system is mathematically unstable.
  • Verify full controllability to confirm a controller can stabilize it.
Derive the linearized state-space matrices

Your nonlinear dynamics function works for simulation, but it is too complex for control design. Linearization computes the Jacobian of the dynamics at the upright equilibrium (theta=0, all velocities zero), giving you matrices A and B that describe how small deviations behave near that point.

  • Create a new cell below your phase portrait cell (Cell 6). Add the linearization code by pasting the following into the new cell:
# Cell 7: Linearization
# Linearized state-space matrices at theta=0 (upright equilibrium)
A = np.array([
    [0, 1, 0, 0],
    [0, 0, -m * g / M, 0],
    [0, 0, 0, 1],
    [0, 0, (M + m) * g / (M * L), 0]
])

B = np.array([
    [0],
    [1 / M],
    [0],
    [-1 / (M * L)]
])

print("Linearized System Matrices:")
print(f"\nA = \n{A}")
print(f"\nB = \n{B.flatten()}")
print(f"\nPhysical meaning:")
print(f"  A[1,2] = {A[1,2]:.4f}: cart acceleration per unit angle (gravity coupling)")
print(f"  A[3,2] = {A[3,2]:.4f}: angular acceleration per unit angle (instability!)")

What does this code do?

The A matrix (4x4) captures how each state variable influences the others near the equilibrium. Each row represents the rate of change of one state.

  • A[1,2] = -m*g/M tells you how much the pendulum angle accelerates the cart (gravity coupling).
  • A[3,2] = (M+m)*g/(M*L) is the key instability term. A positive value here means angle deviations grow exponentially.
  • The B matrix (4x1) captures how the control force enters the system. It affects cart acceleration directly and angular acceleration inversely.
  • Run Cell 7 by pressing Shift+Enter.
  • You should see the A matrix printed as a 4x4 array and B as a 4-element vector, followed by the physical meaning of the key entries.

Seeing a NameError for M, m, L, or g?

The linearization uses the same system parameters defined in Cell 2. Make sure you ran Cell 2 earlier in this session. If your kernel restarted, re-run cells 1 through 6 first.

Why is my linearization cell throwing an error?

Compute eigenvalues to prove instability

The eigenvalues of A are the system's natural frequencies. A positive real eigenvalue means the system diverges exponentially at that rate. This is the mathematical proof of what you saw in the simulation.

  • Create a new cell below Cell 7. Compute and print the eigenvalues by pasting this code:
# Cell 8: Eigenvalue Analysis and Controllability
eigenvalues = np.linalg.eigvals(A)
print("Eigenvalues of A:")
for ev in eigenvalues:
    stability = "UNSTABLE" if ev.real > 0 else "stable" if ev.real < 0 else "marginal"
    print(f"  {ev.real:+.4f} ({stability})")

print(f"\nPositive eigenvalue means the system diverges at rate e^({max(eigenvalues.real):.2f}*t)")
print(f"Time to double: {np.log(2)/max(eigenvalues.real):.3f} seconds")

What do the eigenvalues tell you?

Each eigenvalue represents a natural mode of the linearized system. A positive real part means that mode grows exponentially. A negative real part means it decays.

The "time to double" tells you how fast the instability acts. For this cart-pole, the pendulum's deviation doubles roughly every 0.15 seconds. That is why it crashes so quickly in your earlier simulation.

  • Run Cell 8 by pressing Shift+Enter.
  • You should see four eigenvalues printed, with one labeled UNSTABLE (approximately +4.65). The two zero eigenvalues correspond to the cart's free translation.
Verify controllability

You have proven the system is unstable. The next question is whether a single force input on the cart is enough to control all four states. The controllability matrix answers this definitively.

  • In the same Cell 8, add the controllability check below your eigenvalue code:
# Controllability matrix
n = A.shape[0]
C_ctrl = np.hstack([np.linalg.matrix_power(A, i) @ B for i in range(n)])
rank = np.linalg.matrix_rank(C_ctrl)
print(f"\nControllability matrix rank: {rank} (need {n} for full controllability)")
print(f"System is {'CONTROLLABLE' if rank == n else 'NOT controllable'}!")

What is the controllability matrix?

The controllability matrix C = [B, AB, A^2B, A^3B] tests whether the input can influence all states. If its rank equals the number of states (4), the system is fully controllable.

Full rank means a feedback controller exists that can stabilize the upright equilibrium. The problem HAS a solution. The next step designs that solution.

Before you run this, think: given that the system has one unstable mode and you can only push the cart horizontally, do you think one input is enough to control all four state variables?

  • Run Cell 8 again by pressing Shift+Enter.

You should see the eigenvalue output followed by a controllability matrix rank of 4 and the message "CONTROLLABLE". One positive eigenvalue (approximately +4.65) confirms the instability you observed. Full rank confirms that a feedback controller can fix it.

Seeing rank less than 4?

Double check that your A and B matrices in Cell 7 match the code exactly. A typo in a single entry (especially a sign error) can break the controllability result.

Make sure Cell 7 was run successfully before Cell 8. The controllability code depends on the A and B variables from Cell 7.

Why is my controllability rank not 4?

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with Cells 7 and 8 producing output.

ⓧ I'd like to double check the full code

Your notebook should now have 8 cells in this order. Scroll through and confirm each is present:

  • Cell 1: Imports and Setup (numpy, scipy, matplotlib imports plus print("Setup complete!")).
  • Cell 2: System Parameters (M, m, L, g).
  • Cell 3: Nonlinear Equations of Motion (cart_pole_dynamics function).
  • Cell 4: Simulate Uncontrolled System (sol_uncontrolled).
  • Cell 5: Plot Uncontrolled Response (2x1 subplot, saves uncontrolled_response.png).
  • Cell 6: Phase Portrait (multiple trajectories + quiver plot, saves phase_portrait.png).
  • Cell 7: Linearization (A and B matrices with physical meaning printed).
  • Cell 8: Eigenvalue Analysis and Controllability (eigenvalues with stability labels, time to double, controllability rank).

You now have mathematical proof of two things: the system IS unstable (positive eigenvalue) and it IS controllable (full rank). The problem has a solution. Next up, you will design the LQR controller that actually computes the optimal feedback gains and stabilizes the pendulum.

Design LQR and Stabilize

You have mathematically proven that the pendulum is unstable (positive eigenvalue) and that the system is controllable (full-rank controllability matrix). The pieces are in place for the fix.

Now you will design the controller that actually stabilizes the pendulum. LQR (Linear Quadratic Regulator) finds the optimal feedback gain by solving the continuous-time algebraic Riccati equation. You will close the loop and re-run the same simulation from Step 2. This time, the pendulum stays up.

In this step, get ready to:
  • Design an LQR controller by choosing cost matrices and solving the Riccati equation.
  • Simulate the controlled nonlinear system from the same initial tilt.
  • Create a comparison plot proving the controller stabilizes the pendulum.
Solve the Riccati equation and compute the optimal gain

LQR works by balancing two competing goals: keeping the state small (the pendulum upright, the cart centered) and keeping the control effort small (not slamming the cart back and forth). The Q matrix penalizes state deviations. The R matrix penalizes control effort.

  • Create a new cell in your notebook (Cell 9) and enter this code to design the LQR controller:
# Cell 9: LQR Design
# State cost: penalize angle/angular velocity more than position
Q = np.diag([1.0, 1.0, 10.0, 10.0])
# Control cost
R = np.array([[1.0]])

# Solve the continuous-time algebraic Riccati equation
P = solve_continuous_are(A, B, Q, R)

# Compute optimal feedback gain
K = np.linalg.inv(R) @ B.T @ P
print(f"LQR Gain K = {K.flatten()}")
print(f"\nInterpretation:")
print(f"  K[0] = {K[0,0]:.2f}: gain on cart position")
print(f"  K[1] = {K[0,1]:.2f}: gain on cart velocity")
print(f"  K[2] = {K[0,2]:.2f}: gain on pendulum angle (largest!)")
print(f"  K[3] = {K[0,3]:.2f}: gain on angular velocity")

What does this code do?

Q = diag(1, 1, 10, 10) means you care 10x more about the pendulum angle and angular velocity (indices 2 and 3) than about cart position and cart velocity (indices 0 and 1). This makes physical sense. A wobbly pendulum is dangerous. A slightly off-center cart is fine.

R = 1 sets the price of control effort. A larger R would make the controller gentler but slower. A smaller R would make it aggressive but responsive.

The function solve_continuous_are solves the Riccati equation to find P, from which the optimal gain K is computed as inv(R) @ B.T @ P.

  • Run the cell.
  • You should see the four gain values printed, with K[2] (the angle gain) being the largest.

The gain K tells the controller how strongly to react to each part of the state. A large K[2] means small angular deviations trigger a big corrective push.

But you still need to confirm the closed loop is actually stable. Add this code to the same cell, directly below what you just wrote.

  • Add the closed-loop verification code below the gain printout in Cell 9:
# Verify closed-loop stability
A_cl = A - B @ K
cl_eigenvalues = np.linalg.eigvals(A_cl)
print(f"\nClosed-loop eigenvalues:")
for ev in cl_eigenvalues:
    print(f"  {ev.real:+.4f} + {ev.imag:+.4f}j")
print("All negative real parts = STABLE!")

Why check the closed-loop eigenvalues?

In Step 4, you saw the open-loop system had a positive eigenvalue (~+4.65), proving instability. Now A_cl = A - B @ K is the closed-loop system matrix. If all its eigenvalues have negative real parts, the feedback has moved the system from unstable to stable.

This is the mathematical proof that your controller works, before you even simulate it.

  • Re-run Cell 9 (the full cell now includes both the gain computation and the verification).
  • You should see four closed-loop eigenvalues, all with negative real parts, followed by the confirmation message.

Seeing an error about A or B not defined?

Make sure Cell 7 (where you defined the linearized matrices A and B) has been executed in this notebook session. If you restarted the kernel, re-run all cells from Cell 1 through Cell 8 first.

Help me figure out which prior cells I need to re-run.

Simulate the controlled nonlinear system

You have the gain K and mathematical proof it works on the linearized model. Now you apply it to the full nonlinear dynamics from Step 2 to see if it works on the real system.

  • Create Cell 10 and enter this code to simulate the controlled system:
# Cell 10: Simulate Controlled System
def controlled_dynamics(t, state):
    """Nonlinear dynamics with LQR feedback."""
    u = float(-K @ np.array(state).reshape(-1, 1))
    u = np.clip(u, -50, 50)  # realistic force limit (Newtons)
    return cart_pole_dynamics(t, state, F=u)

sol_controlled = solve_ivp(
    controlled_dynamics, t_span, state0,
    t_eval=t_eval, method='RK45'
)

# Compute control effort for plotting
control_effort = np.array([
    float(-K @ sol_controlled.y[:, i].reshape(-1, 1))
    for i in range(len(sol_controlled.t))
])
control_effort = np.clip(control_effort, -50, 50)

What does this code do?

The controlled_dynamics function wraps your existing cart_pole_dynamics but computes u = -K @ state at every timestep. This is full state feedback: the controller measures the entire state (position, velocity, angle, angular velocity) and pushes the cart to cancel any deviation.

The np.clip(u, -50, 50) enforces a realistic actuator limit of 50 Newtons. A real motor cannot apply infinite force.

The simulation uses the same t_span, state0, and t_eval from Step 2, so you can compare directly.

  • Run Cell 10.
  • The cell completes without error. No plot appears yet. The next cell will visualize the results.

Seeing a ValueError about array shapes?

Make sure K is a 1x4 matrix (from Cell 9) and that state is being reshaped to a column vector with .reshape(-1, 1). If K has the wrong shape, re-run Cell 9.

Help me debug the matrix dimensions.

Create the comparison figure

This is the payoff. You will overlay the controlled and uncontrolled trajectories on the same axes, proving the difference feedback makes.

  • Create Cell 11 and enter this code to build the comparison plot:
# Cell 11: Comparison Plot
fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)

# Top: angle comparison (uncontrolled vs controlled)
axes[0].plot(sol_uncontrolled.t, np.degrees(sol_uncontrolled.y[2]),
             'r--', linewidth=2, label='Uncontrolled')
axes[0].plot(sol_controlled.t, np.degrees(sol_controlled.y[2]),
             'b', linewidth=2, label='LQR Controlled')
axes[0].axhline(y=0, color='k', linestyle=':', alpha=0.3)
axes[0].set_ylabel('Pendulum Angle (degrees)')
axes[0].set_title('LQR Stabilization of the Inverted Pendulum')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Middle: cart position under control
axes[1].plot(sol_controlled.t, sol_controlled.y[0], 'b', linewidth=2)
axes[1].set_ylabel('Cart Position (m)')
axes[1].grid(True, alpha=0.3)

# Bottom: control force applied
axes[2].plot(sol_controlled.t, control_effort, 'g', linewidth=2)
axes[2].set_xlabel('Time (s)')
axes[2].set_ylabel('Control Force (N)')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lqr_stabilization.png', dpi=150)
plt.show()

What does this plot show?

The top subplot overlays the uncontrolled angle (red dashed, diverging past 90 degrees) against the controlled angle (blue solid, converging back to zero). This is the core proof that LQR works.

The middle subplot shows the cart moving slightly to push the pendulum back. The bottom subplot reveals the control force required. You will see a spike at the start as the controller catches the tilt, then it settles to near zero.

Before you run this cell, predict: how quickly do you think the controlled angle will return to zero? Will it overshoot past zero before settling?

  • Run Cell 11.
  • You should see a 3-subplot figure. In the top subplot, the blue controlled line converges to zero within about 2 seconds while the red uncontrolled line diverges wildly. The middle subplot shows the cart settling. The bottom subplot shows the control force.

The pendulum stays upright. This is the same system that crashed within 1.5 seconds in Step 2, now stabilized indefinitely by feedback.

Plot looks empty or the controlled line diverges?

If the controlled line also diverges, check that Cell 9 ran successfully and that K has sensible values (K[2] should be around -40 to -60). Re-run Cells 9 and 10 before re-running Cell 11.

If the plot is blank, make sure sol_controlled and sol_uncontrolled both exist by re-running Cells 4 and 10.

Help me debug my LQR comparison plot.

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook with all cells from 1 through 11 executed.

ⓧ I'd like to double check the full code

Your notebook should now contain 11 cells. Scroll through and confirm these three new cells are present in order:

  • Cell 9 (LQR Design): defines Q, R, solves the Riccati equation, computes K, prints the gain interpretation, and verifies closed-loop stability.
  • Cell 10 (Simulate Controlled System): defines controlled_dynamics, runs solve_ivp, and computes control_effort.
  • Cell 11 (Comparison Plot): creates the 3-subplot figure and saves lqr_stabilization.png.

You have gone from a system that crashes in 1.5 seconds to one that balances indefinitely. In the Secret Mission, you will push this controller harder by throwing random disturbances at the pendulum and tuning Q/R to explore the trade-off between fast recovery and smooth control.

Secret mission

Disturbance Rejection and Q/R Tuning

Your LQR controller stabilizes the pendulum from a gentle initial tilt. But a real robot gets pushed, bumped, and kicked. Can your controller handle random impulse disturbances mid-flight? And what happens when you crank up the aggressiveness of the control response?

Clean Up Your Resources

Clean Up Your Resources

Decide whether to keep your resources running, pause them to come back later, or delete them entirely. This project runs entirely locally, so there are no ongoing costs.

Resources you used:

  • The control-lab virtual environment (contains NumPy, SciPy, Matplotlib, and Jupyter Notebook).
  • The cart_pole_control.ipynb notebook file.
  • The cart_pole_control_secret.py secret mission file.
  • Generated plot files: uncontrolled_response.png, phase_portrait.png, lqr_stabilization.png, and disturbance_rejection.png.

Keep everything running

No action needed. Choose this if you plan to keep experimenting with control theory, tuning Q/R matrices, or extending the notebook with new controllers.

  • Your control-lab virtual environment stays on disk and you can reactivate it any time.
  • Your notebook and generated plots remain in place for future reference.

Pause - I'll come back to this later

Shut down Jupyter and deactivate the virtual environment to free up memory. Your files stay on disk so you can pick up where you left off.

  • In the terminal where Jupyter is running, press Ctrl+C to shut down the notebook server.
  • Deactivate the virtual environment by running:
deactivate
  • To resume later, navigate back to your Desktop, run source control-lab/bin/activate, and launch jupyter notebook again.

Delete - I don't want to use this again

Remove all project resources and start fresh. This deletes the virtual environment, notebook, script, and generated plots.

  • In the terminal where Jupyter is running, press Ctrl+C to shut down the notebook server.
  • Deactivate the virtual environment by running:
deactivate
  • Navigate to your Desktop and remove the virtual environment by running:
cd ~/Desktop
rm -rf control-lab
  • Remove the notebook, script, and generated plot files by running:
rm cart_pole_control.ipynb cart_pole_control_secret.py
rm uncontrolled_response.png phase_portrait.png lqr_stabilization.png disturbance_rejection.png
  • Confirm everything is removed by running:
ls
  • You should no longer see control-lab, cart_pole_control.ipynb, or any of the .png files in the output.

That's a Wrap!

That's a Wrap!

Nice work. You just built a complete control design pipeline for one of the most iconic problems in robotics: stabilizing an inverted pendulum.

You've learned how to:

  • Simulate the full nonlinear dynamics of a cart-pole system from first principles using the Euler-Lagrange equations and SciPy's solve_ivp.
  • Visualize instability through phase portraits and prove it mathematically with eigenvalue analysis and the controllability rank test.
  • Design an LQR controller by solving the Riccati equation and demonstrate stabilization of the unstable upright equilibrium.
  • (Secret Mission) Implemented disturbance rejection with random impulse kicks and explored the Q/R tuning trade-off that governs every real-world controller.

Ready to quiz yourself?