🧠 Continuous-Time Luenberger Observer Design in Python

In this tutorial, we demonstrate how to design a Luenberger observer for a continuous-time system using Python. This example uses a DC motor modeled in state space and simulates how an observer estimates internal system states from outputs.

🔧 Step 1: Define the System

A DC motor is a mechatronic system containing both electrical and mechanical components. It can be modeled as the following:

\(\frac{d}{dt} \begin{bmatrix} i \\ \theta \\ \dot{\theta} \end{bmatrix} = \begin{bmatrix} -\frac{R}{L} & 0 & -\frac{K}{L} \\ 0 & 0 & 1 \\ \frac{K}{J} & 0 & -\frac{B}{J} \end{bmatrix} \begin{bmatrix} i \\ \theta \\ \dot{\theta} \end{bmatrix} + \begin{bmatrix} \frac{1}{L} \\ 0 \\ 0 \end{bmatrix} V \)

\( y = \begin{bmatrix} 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} i \\ \theta \\ \dot{\theta} \end{bmatrix} \)

🟢 Stability Analysis: Eigenvalues of A

Let us first take a look at the system stability by checking the eigenvalues of the “A” matrix:

eigvals = np.linalg.eig(A)
print("Original Eigenvalues:", eigvals[0])

Execute the code and see what you will get.

🔍 Step 2: Check Observability

Before designing the observer, we check that the system is observable:

O = np.linalg.matrix_rank(np.concatenate((C, C @ A, C @ A @ A), axis=0))
print("Observability Rank:", O)

You can also use the np.block function to construct the observability matrix:

O = np.linalg.matrix_rank(np.block([ [C], [C@A], [C@(A@A)] ]))
print("Observability Rank:", O)

You should see the output to be similar to the following:

Observability Rank: 3

The full rank of the observability matrix indicates that the system is observable.

🏀 Step 3: Design the Observer

To design the observer, we select the following desired poles for the estimator error dynamics:

import control.matlab as ct

pole_des = np.array([-500+250j, -500-250j, -1000])
L = ct.place(A.T, C.T, pole_des).T

This way, the matrix \( L \) will ensure that the observer dynamics \( A - LC \) have the desired fast convergence behavior. Let us verify that this is indeed the case by computing the actual eigenvalues of \( A - LC \):

est_poles = np.linalg.eig(A - L @ C)
print("Observer Poles:", est_poles[0])

You should see the following:

Observer Poles: [-500.+250.j -500.-250.j -1000.+0.j]

♻️ Step 4: Simulate the Observer

To simulate the observer, we create an augmented system that includes both the plant and the observer:

# Augmented system
Aaug = np.block([[A, np.zeros((3, 3))], [L @ C, A - L @ C]])
Baug = np.vstack((B, B))
Caug = np.hstack((C, np.zeros((1, 3))))
Daug = np.array([0])
sys = ct.ss(Aaug, Baug, Caug, Daug)


x0 = np.array([10, 2, 10]); xhat0 = np.array([0, 0, 0]); X0 = np.array([x0, xhat0]).reshape((6, 1))

Tend = 0.03; amplitude = 10; initpha = 0; freq = 600
t = np.arange(0, Tend, 1e-4)

u = amplitude * np.sin(freq * t + initpha)

[Y, T, X] = ct.lsim(sys, u, t, X0)

📊 Step 5: Plot the Results

import matplotlib.pyplot as plt

plt.figure()
plt.subplot(3, 1, 1)
plt.plot(t, X[:, 0], t, X[:, 3], '--', linewidth=1.5)
plt.xlabel('time (sec)')
plt.legend(['$x_1 = i_a$', '$\hat x_1$'], fontsize=16)
plt.grid()
plt.ylabel('$x_1$', fontsize=16)
plt.title('States and their estimates')
plt.subplot(3, 1, 2)
plt.plot(t, X[:, 1], t, X[:, 4], '--', linewidth=1.5)
plt.xlabel('time (sec)')
plt.legend(['$x_2 = \\theta$', '$\hat x_2$'], fontsize=16)
plt.grid()
plt.ylabel('$x_2$', fontsize=16)
plt.subplot(3, 1, 3)
plt.plot(t, X[:, 2], t, X[:, 5], '--', linewidth=1.5)
plt.xlabel('time (sec)')
plt.legend(['$x_3 = \dot{\\theta}$', '$\hat x_3$'], fontsize=16)
plt.grid()
plt.ylabel('$x_3$', fontsize=16)
plt.show()

You should see something similar to the following, where the observer accurately tracks the angular position and velocity over time, with rapid convergence.

Last updated on May 9, 2025