# Example for Linear Quadratic Gaussion control design
**URL:** https://discourse.julialang.org/t/example-for-linear-quadratic-gaussion-control-design/11829
**Category:** General Usage
**Tags:** control
**Created:** [June 20, 2018, 3:43pm UTC](https://discourse.julialang.org/t/example-for-linear-quadratic-gaussion-control-design/11829 "2018-06-20T15:43:34Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [June 20, 2018, 3:43pm UTC](https://discourse.julialang.org/t/example-for-linear-quadratic-gaussion-control-design/11829/1 "2018-06-20T15:43:34Z")
Hello,
My colleague and me wrote the following example for LQG design of a speed controller of a DC motor:
```julia
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.
```julia
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
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [December 4, 2022, 9:19am UTC](https://discourse.julialang.org/t/example-for-linear-quadratic-gaussion-control-design/11829/2 "2022-12-04T09:19:48Z")
Here’s how your example would look now
```julia
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`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/api/#RobustAndOptimalControl.closedloop) for more details. See also [`ExtendedStateSpace`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/api/#RobustAndOptimalControl.ExtendedStateSpace) for how the two-input, two-output system depicted below works
```julia
z ┌─────┐ w
◄──┤ │◄──
│ P │
┌──┤ │◄─┐
│y └─────┘ u│
│ │
│ ┌─────┐ │
│ │ │ │
└─►│ K ├──┘
│ │
└─────┘
```
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [December 4, 2022, 6:24pm UTC](https://discourse.julialang.org/t/example-for-linear-quadratic-gaussion-control-design/11829/3 "2022-12-04T18:24:21Z")
Thank you very much! 😀