Example for Linear Quadratic Gaussion control design

ufechner7 · June 20, 2018, 3:43pm

Hello,

My colleague and me wrote the following example for LQG design of a speed controller of a DC motor:

using ControlSystems
"""
Example for designing an LQG speed controller for an electrical DC motor.
"""

# Constants
Ke = 0.006156                                # electromotive force constant in V/rpm
Kt                      = 0.0728             # Torque constant (Nm/A)
J                       = 2.8*700e-6;        # Inertia of motor plus load (kgm^2)
Rel                     = 0.11;              # DC motor resistance (Ω)
L                       = 34e-6;             # DC motor inductance (H)

# helper to convert a scalar into an 1x1 Matrix
mat(sc) = reshape([sc],1,1)

# Create an instance of a DC motor model in state space form.
# Ke: electromotive force constant in V/rpm
# Kt: torque constant (Nm/A)
# L: Inductance; R: Resistance, J: Inertia; b: viscous friction coefficient
function motor(Ke, Kt, L, R, J, b=1e-3)
    Ke1 = Ke*60.0/2π   # change unit to rad/s
    A = [-b/J     Kt/J
         -Ke1/L   -R/L];
    B = [0
        1/L];
    C = [1     0];
    D = 0;
    sys = ss(A, B, C, D);
end

p60 = motor(Ke, Kt, L, Rel, J)
stepplot(p60, 0.2, 0.001)
bodeplot(p60)

# LQR control
Q = [1.     0;
     0      1]

R = 20.
K = lqr(p60.A, p60.B, Q, R)
# needs to be modified if Nbar is not a scalar
Nbar = 1. / (p60.D - (p60.C - p60.D*K) * inv(p60.A - p60.B*K) * p60.B)

# Kalman filter based observer
Vd = [10.     0   # covariance of the speed estimation
     0       10]; # covariance of the current estimation
Vn = 0.04;       # covariance for the speed measurement (radians/s)^2
G = LQG(p60, Q, mat(R), Vd, mat(Vn))
Gcl = G[:cl]
T = G[:T]
S = G[:S];
f1 = sigmaplot([S,T],logspace(-3,3,1000))
f2 = stepplot(Gcl, label=["Closed loop system using LQG"])
Plots.plot(f1,f2)

It requires to use the master branch of ControlSystems.jl.

Pkg.add("ControlSystems.jl")
Pkg.checkout("ControlSystems.jl")

Any suggestions how to simplify or improve this code?
For example, is the gain Nbar also available using the LQG type?

Extra question: Is it possible to get more then one plot window?

Uwe

baggepinnen · December 4, 2022, 9:19am

Here’s how your example would look now

using ControlSystemsBase, RobustAndOptimalControl, LinearAlgebra, Plots

# Constants
Ke = 0.006156                                # electromotive force constant in V/rpm
Kt                      = 0.0728             # Torque constant (Nm/A)
J                       = 2.8*700e-6;        # Inertia of motor plus load (kgm^2)
Rel                     = 0.11;              # DC motor resistance (Ω)
L                       = 34e-6;             # DC motor inductance (H)


# Create an instance of a DC motor model in state space form.
# Ke: electromotive force constant in V/rpm
# Kt: torque constant (Nm/A)
# L: Inductance; R: Resistance, J: Inertia; b: viscous friction coefficient
function motor(Ke, Kt, L, R, J, b=1e-3)
    Ke1 = Ke*60.0/2π   # change unit to rad/s
    A = [-b/J     Kt/J
         -Ke1/L   -R/L];
    B = [0
        1/L];
    C = [1     0];
    D = 0;
    sys = ss(A, B, C, D);
end

p60 = motor(Ke, Kt, L, Rel, J)
plot(step(p60, 0:0.001:0.2, method=:zoh))
bodeplot(p60)

# LQR control
Q = [1.     0;
     0      1]

R = 20.0I(1)
K = lqr(p60, Q, R)
# needs to be modified if N̄ is not a scalar
N̄ = [1.] / (p60.D - (p60.C - p60.D*K) * inv(p60.A - p60.B*K) * p60.B)

# Kalman filter based observer
Vd = [10.     0   # covariance of the speed estimation
     0       10]; # covariance of the current estimation
Vn = 0.04I(1);       # covariance for the speed measurement (radians/s)^2
# G = LQGProblem(ExtendedStateSpace(p60, C1=I(2), B1=I(2)), Q, R, Vd, Vn)
G = LQGProblem(ExtendedStateSpace(p60, C1=I(2), B1=I(2)), Q, R, Vd, Vn)
Gclwz = lft(G) # Transfer function from disturbance inputs to controlled outputs
Gcl = closedloop(G) # Transfer function from reference inputs to controlled outputs

T = output_comp_sensitivity(G) 
S = output_sensitivity(G);

logspace(l,u,n) = exp10.(range(l, u, length=n))
f1 = sigmaplot([S,T],logspace(-3,4,1000))
f2 = plot(step(T*inv(dcgain(T)[]), 0:0.01:2, method=:zoh), title=["Closed loop system using LQG"])
f2 = plot(step(Gcl*N̄[], 0:0.01:2, method=:zoh), title=["Closed loop system using LQG"])
Plots.plot(f1,f2)

using Test
@test static_gain_compensation(ssdata(p60)..., lqr(G)) ≈ N̄
dcgain(Gcl*static_gain_compensation(G))

Different authors define the “closed-loop transfer function” in the LQG context differently, the function static_gain_compensation follows the definition of Gland and Ljung, see closedloop for more details. See also ExtendedStateSpace for how the two-input, two-output system depicted below works

z  ┌─────┐  w
◄──┤     │◄──
   │  P  │
┌──┤     │◄─┐
│y └─────┘ u│
│           │
│  ┌─────┐  │
│  │     │  │
└─►│  K  ├──┘
   │     │
   └─────┘

ufechner7 · December 4, 2022, 6:24pm

Thank you very much!