Tutorial: Control

Tutorial: Control#

Control I – Classical Feedback#

To support reusability in other projects, here’s a clean, general PID controller model implemented as a Python class. This can be copied into notebooks for simulating PID in various systems (e.g., orbital control or thermal regulation). It supports discrete-time approximation for numerical integration, making it versatile for custom simulations.

import numpy as np

class PIDController:
    """General PID controller class for simulation.

    Parameters:
    - Kp: Proportional gain
    - Ki: Integral gain
    - Kd: Derivative gain
    - dt: Time step for discrete approximation (default 0.01 s)

    Usage:
    pid = PIDController(Kp=1.0, Ki=0.5, Kd=0.1, dt=0.01)
    control_signal = pid.update(setpoint, measurement)
    """
    def __init__(self, Kp=0.0, Ki=0.0, Kd=0.0, dt=0.01):
        self.Kp = Kp
        self.Ki = Ki
        self.Kd = Kd
        self.dt = dt
        self.integral = 0.0
        self.prev_error = 0.0

    def update(self, setpoint, measurement):
        error = setpoint - measurement
        self.integral += error * self.dt
        derivative = (error - self.prev_error) / self.dt
        output = self.Kp * error + self.Ki * self.integral + self.Kd * derivative
        self.prev_error = error
        return output

    def reset(self):
        self.integral = 0.0
        self.prev_error = 0.0

Exercise 1: Laplace Transform for Spacecraft Attitude Dynamics (30-45 minutes)#

Use Laplace transforms to solve an ODE for satellite attitude and verify numerically.

Model satellite rotation as \(J \frac{d^2 \theta}{dt^2} + c \frac{d\theta}{dt} = \tau(t)\) (J=10 kg·m² inertia, c=0.5 N·m·s/rad damping, \(\tau(t)\)=step torque of 1 N·m).

  1. Derive the transfer function \(\Theta(s)/T(s)\) symbolically with sympy.

  2. Compute inverse Laplace for \(\theta(t)\) and plot analytically.

  3. Simulate with control.step_response and compare curves.

  4. Discuss: How does damping affect stability in real satellite pointing?

Starter Code:

import sympy as sp
import control as ct
import matplotlib.pyplot as plt
import numpy as np

s, t = sp.symbols('s t')
J, c = 10, 0.5
T_s = 1 / s  # Step input

# Derive Theta(s)...

theta_t = sp.inverse_laplace_transform(..., s, t)
print(theta_t)

G = ct.tf([1], [J, c, 0])
t_sim = np.linspace(0, 50, 1000)
_, theta_sim = ct.step_response(G, t_sim)

# Plot

---------------------------------------------------------------------------
SympifyError                              Traceback (most recent call last)
Cell In[2], line 12
      8 T_s = 1 / s  # Step input
     10 # Derive Theta(s)...
---> 12 theta_t = sp.inverse_laplace_transform(..., s, t)
     13 print(theta_t)
     15 G = ct.tf([1], [J, c, 0])

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/integrals/laplace.py:2326, in inverse_laplace_transform(F, s, t, plane, **hints)
   2322 if isinstance(F, MatrixBase) and hasattr(F, 'applyfunc'):
   2323     return F.applyfunc(
   2324         lambda Fij: inverse_laplace_transform(Fij, s, t, plane, **hints))
-> 2326 r, c = InverseLaplaceTransform(F, s, t, plane).doit(
   2327     noconds=False, simplify=_simplify)
   2329 if _noconds:
   2330     return r

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/integrals/laplace.py:2220, in InverseLaplaceTransform.__new__(cls, F, s, x, plane, **opts)
   2218 if plane is None:
   2219     plane = InverseLaplaceTransform._none_sentinel
-> 2220 return IntegralTransform.__new__(cls, F, s, x, plane, **opts)

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
     69 @wraps(func)
     70 def wrapper(*args, **kwargs):
     71     try:
---> 72         retval = cfunc(*args, **kwargs)
     73     except TypeError as e:
     74         if not e.args or not e.args[0].startswith('unhashable type:'):

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/core/function.py:450, in Function.__new__(cls, *args, **options)
    448     return UndefinedFunction(*args, **options)  # type: ignore
    449 else:
--> 450     return cls._new_(*args, **options)

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/core/function.py:472, in Function._new_(cls, *args, **options)
    464     raise TypeError(temp % {
    465         'name': cls,
    466         'qual': 'exactly' if len(cls.nargs) == 1 else 'at least',
    467         'args': min(cls.nargs),
    468         'plural': 's'*(min(cls.nargs) != 1),
    469         'given': n})
    471 evaluate = options.get('evaluate', global_parameters.evaluate)
--> 472 result = super().__new__(cls, *args, **options)
    473 if evaluate and isinstance(result, cls) and result.args:
    474     _should_evalf = [cls._should_evalf(a) for a in result.args]

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
     69 @wraps(func)
     70 def wrapper(*args, **kwargs):
     71     try:
---> 72         retval = cfunc(*args, **kwargs)
     73     except TypeError as e:
     74         if not e.args or not e.args[0].startswith('unhashable type:'):

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/core/function.py:299, in Application.__new__(cls, *args, **options)
    296 from sympy.sets.fancysets import Naturals0
    297 from sympy.sets.sets import FiniteSet
--> 299 args = list(map(sympify, args))
    300 evaluate = options.pop('evaluate', global_parameters.evaluate)
    301 # WildFunction (and anything else like it) may have nargs defined
    302 # and we throw that value away here

File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/sympy/core/sympify.py:491, in sympify(a, locals, convert_xor, strict, rational, evaluate)
    488         pass
    490 if not isinstance(a, str):
--> 491     raise SympifyError('cannot sympify object of type %r' % type(a))
    493 from sympy.parsing.sympy_parser import (parse_expr, TokenError,
    494                                         standard_transformations)
    495 from sympy.parsing.sympy_parser import convert_xor as t_convert_xor

SympifyError: SympifyError: "cannot sympify object of type <class 'ellipsis'>"

Exercise 2: Control Loop Anatomy and Basic Feedback Simulation (30-45 minutes)#

Build and simulate a feedback loop for orbital adjustment.

  1. For plant \(G(s) = \frac{1}{s^2 + 0.5s + 1}\) (relative orbit with drag), simulate open/closed-loop (Kp=5) step responses.

  2. Add disturbance (e.g., constant bias in input) and recompute.

  3. Use control.step_info for metrics.

  4. Discuss: Feedback’s role in rejecting space disturbances like solar wind.

Starter Code:

import control as ct
import matplotlib.pyplot as plt
import numpy as np

G = ct.tf([1], [1, 0.5, 1])
t = np.linspace(0, 20, 1000)

# Open-loop...

# Closed-loop...

# Plot and metrics...

Exercise 3: Time-Domain Simulation vs. Control Package Comparison (45-60 minutes)#

Implement manual time-domain simulation and compare to control package for PID validation.

  1. Use the PIDController class above to simulate a mass-spring-damper plant (\(m=1\), \(c=0.5\), \(k=1\)) with PID (Kp=5, Ki=1, Kd=1) via Euler integration (discrete steps).

  2. Compare to control.step_response output for the same system.

  3. Plot both responses and compute differences (e.g., RMSE).

  4. Add noise to manual simulation and discuss numerical accuracy in space control (e.g., onboard computation limits).

Starter Code:

import numpy as np
import control as ct
import matplotlib.pyplot as plt

# PID class from above...

dt = 0.01
t = np.arange(0, 20, dt)
setpoint = np.ones_like(t)  # Step

# Manual simulation (Euler)
y = 0; dy = 0
pid = PIDController(Kp=5, Ki=1, Kd=1, dt=dt)
ys_man = []
for sp in setpoint:
    u = pid.update(sp, y)
    ddy = (u - 0.5 * dy - y) / 1  # Plant ODE
    dy += ddy * dt
    y += dy * dt
    ys_man.append(y)

# Control package
G = ct.tf([1], [1, 0.5, 1])
C = ct.tf([1, 5, 1], [1, 0])  # Kd s^2 + Kp s + Ki / s ? Wait, adjust
sys = ct.feedback(C * G, 1)
_, ys_ctrl = ct.step_response(sys, t)

# Plot comparison, RMSE = np.sqrt(np.mean((ys_man - ys_ctrl)**2))

Exercise 4: Stability Margins and Bode Tuning (30-45 minutes)#

Tune using Bode plots for robustness in flexible structures.

  1. For plant \(G(s) = \frac{10}{s(s+2)(s+5)}\), plot Bode for varying Kp and compute margins.

  2. Tune for PM > 45°.

  3. Discuss: Margins’ importance in satellite flexible mode suppression.

Starter Code:

import control as ct
import matplotlib.pyplot as plt

G = ct.tf([10], [1, 7, 10, 0])

# Bode for Kp=...

gm, pm, _, _ = ct.margin(Kp * G)

Exercise 5: Ziegler-Nichols Application and Variants (30-45 minutes)#

Apply ZN and compare variants.

  1. Find Ku, Pu for the plant in Exercise 4.

  2. Compute gains for all variants and plot responses.

  3. Discuss: Variant choice for low-overshoot space applications.

Starter Code:

# From chapter ZN code...

Control II – State-Space & Optimal Linear Control#

Excercise 1: State-Space Representation and System Properties#

Consider a simplified satellite attitude dynamics model (double integrator) from Section 11.2: second-order system with transfer function \( G(s) = \frac{1}{I s^2} \) (moment of inertia \( I=1 \)), input torque \( u = \tau \), states angle \( \theta = y_1 \) and rate \( \dot{\theta} = y_2 \), output angle (\( \mathbf{C} = [1, 0] \), \( \mathbf{D} = 0 \)).

  1. Derive the state-space matrices \( \mathbf{A} \), \( \mathbf{B} \).

  2. Compute controllability and observability matrices using NumPy; check ranks and singular values.

  3. Simulate the step response (\( u=1 \)) for 10 seconds using solve_ivp, plotting states. Add a small damping term to \( \mathbf{A} \) (e.g., -0.1 on \( A_{22} \)) and re-analyze stability via eigenvalues.

  4. Discuss: How does uncontrollability arise in a faulty thruster scenario?

Hint: Use the plot_step_response function from the chapter for visualization.

Excercise 2: MIMO Decoupling and RGA Analysis (20 minutes)#

Use the 2x2 coupled system from Section 11.3.5: \( \mathbf{A} = \begin{bmatrix} -1 & 0.5 \\ 0.5 & -2 \end{bmatrix} \), \( \mathbf{B} = \begin{bmatrix} 1 & 0.2 \\ 0.1 & 1 \end{bmatrix} \), \( \mathbf{C} = \mathbf{I}_2 \), \( \mathbf{D} = \mathbf{0} \).

  1. Compute steady-state gain \( \mathbf{G} = -\mathbf{C} \mathbf{A}^{-1} \mathbf{B} \) using NumPy.

  2. Calculate RGA \( \mathbf{\Lambda} = \mathbf{G} \odot (\mathbf{G}^{-1})^\top \); interpret pairings (diagonal >1 indicates interactions).

  3. Design static decoupler \( \mathbf{D}_s = \mathbf{G}^{-1} \); simulate step on first input with/without decoupler using solve_ivp, plotting outputs to show reduced cross-coupling.

  4. Discuss: In a spacecraft with coupled thrusters, how does decoupling improve independent axis control?

Excercise 3: LQR Design and Tuning#

Objective: Design and tune an LQR controller, verify performance. Use the mass-spring-damper model from Section 11.2 (\( \omega_n=1 \), \( \zeta=0.5 \)): \( \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -1 & -1 \end{bmatrix} \), \( \mathbf{B} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \).

(a) Solve the continuous ARE for infinite-horizon LQR with \( \mathbf{Q} = \mathbf{I}_2 \), \( \mathbf{R} = [[0.1]] \); compute \( \mathbf{K} \) using solve_continuous_are. (b) Simulate closed-loop response from initial \( [1, 0]^\top \) for 10 seconds; plot states and control input \( u = -\mathbf{K} \mathbf{y} \). (c) Tune: Increase \( \mathbf{Q}_{11} \) to 10 (prioritize position); re-compute \( \mathbf{K} \) and simulate. Compare settling times and input magnitudes. (d) Verify closed-loop stability with eigenvalues of \( \mathbf{A} - \mathbf{B} \mathbf{K} \). (e) Extend to space: Interpret as vibration control in a flexible spacecraft appendage; discuss fuel implications of high \( \mathbf{R} \).

Hint: Integrate the input computation in the ODE function for solve_ivp.

Problem 3: Kalman Filter Implementation (30 minutes) Objective: Implement and simulate a discrete-time Kalman filter for state estimation. Use the discretized satellite model ( \( \Delta t = 0.1 \) ): approximate \( \mathbf{A}_d = \mathbf{I} + \mathbf{A} \Delta t \), \( \mathbf{B}_d = \mathbf{B} \Delta t \), measuring angle only (\( \mathbf{C} = [1, 0] \)).

(a) Simulate true noisy trajectory for 100 steps from \( [1, 0]^\top \), with \( \mathbf{Q}_{n,d} = 0.01 \mathbf{I}_2 \), \( \mathbf{R}_n = 0.1 \); generate measurements \( z_k \). (b) Implement the Kalman filter loop: prediction and update steps, starting from \( \hat{\mathbf{y}}_0 = [0, 0]^\top \), \( \mathbf{P}_{e,0} = \mathbf{I}_2 \). (c) Plot true states, estimates, and noisy measurements. Compute root-mean-square estimation error. (d) Tune: Increase \( \mathbf{R}_n \) to 1 (noisier sensor); re-simulate and discuss filter reliance on model vs. measurements. (e) Relate to space: How could this estimate angular rate from noisy star-tracker data?

Hint: Use the discrete-time equations from Section 11.5.1; precompute steady-state \( \mathbf{L} \) if time allows.

Excercise 4: LQR Design and Tuning#

Use the mass-spring-damper model from Section 11.2 (\( \omega_n=1 \), \( \zeta=0.5 \)): \( \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -1 & -1 \end{bmatrix} \), \( \mathbf{B} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \).

  1. Solve the continuous ARE for infinite-horizon LQR with \( \mathbf{Q} = \mathbf{I}_2 \), \( \mathbf{R} = [[0.1]] \); compute \( \mathbf{K} \) using solve_continuous_are.

  2. Simulate closed-loop response from initial \( [1, 0]^\top \) for 10 seconds; plot states and control input \( u = -\mathbf{K} \mathbf{y} \).

  3. Tune: Increase \( \mathbf{Q}_{11} \) to 10 (prioritize position); re-compute \( \mathbf{K} \) and simulate. Compare settling times and input magnitudes.

  4. Verify closed-loop stability with eigenvalues of \( \mathbf{A} - \mathbf{B} \mathbf{K} \).

  5. Extend to space: Interpret as vibration control in a flexible spacecraft appendage; discuss fuel implications of high \( \mathbf{R} \).

(Optional) Excercise 5: Integrated LQG for Noisy System#

Build on the mass-spring-damper (exercise 3): Add noise \( \mathbf{Q}_n = 0.01 \mathbf{I}_2 \), \( \mathbf{R}_n = 0.1 \), partial output \( \mathbf{C} = [1, 0] \).

  1. Design LQR (\( \mathbf{Q} = \mathbf{I}_2 \), \( \mathbf{R} = 0.1 \)) and Kalman separately.

  2. Implement augmented LQG dynamics in solve_ivp; simulate from initial offset with noise, plotting true/estimated states.

  3. Verify separation: Compare closed-loop eigenvalues to individual LQR/Kalman ones.

  4. Tune for space scenario: Increase \( \mathbf{R} \) for “fuel-efficient” control; discuss estimation convergence under high noise.

Hint: Use the provided LQG example code as a base; focus on interpretation.

Contents