Nonlinear functionality
The nonlinear interface is currently experimental and at any time subject to breaking changes not respecting semantic versioning.
ControlSystems.jl can represent nonlinear feedback systems that can be written on the form
┌─────────┐
y◄───┤ │◄────u
│ P │
Δy┌───┤ │◄───┐Δu
│ └─────────┘ │
│ │
│ ┌───┐ │
└─────►│ f ├───────┘
└───┘i.e., as a linear-fractional transform (LFT) between a linear system $P$ and a diagonal matrix with scalar non-linear functions $f$. This representation is identical to that used for delay systems, and is exposed to the user in a similar way as well. The main entry point is the function nonlinearity which takes a nonlinear function $f$ like so, nonlinearity(f). This creates a primitive system containing only the nonlinearity, but which behaves like a standard LTISystem during algebraic operations. We illustrate its usage through a number of examples.
Examples
Control-signal saturation
To create a controller that saturates the output at $\pm 0.7$, we call
using ControlSystems, Plots
using ControlSystemsBase: nonlinearity # This functionality is not exported due to the beta status
C = pid(1, 0.1, form=:parallel) # A standard PI controller
nl = nonlinearity(x->clamp(x, -0.7, 0.7)) # a saturating nonlinearity
satC = nl*C # Connect the saturation at the output of CControlSystemsBase.HammersteinWienerSystem{Float64}
P: StateSpace{Continuous, Float64}
A =
0.0
B =
0.25 0.0
C =
0.0
0.4
D =
0.0 1.0
1.0 0.0
Continuous-time state-space model
Nonlinearities: Function[Main.var"#1#2"()]we may now use this controller like we would normally do in ControlSystems, e.g.,
P = tf(1, [1, 1]) # a plant
G = feedback(P*C) # closed loop without nonlinearity
Gnl = feedback(P*satC) # closed loop with saturation
Gu = feedback(C, P) # closed loop from reference to control signal without nonlinearity
Gunl = feedback(satC, P) # closed loop from reference to control signal with saturation
plot(step([G; Gu], 5), lab = ["Linear y" "Linear u"])
plot!(step([Gnl; Gunl], 5), lab = ["Nonlinear y" "Nonlinear u"])Since the saturating nonlinearity is common, we provide the constructor ControlSystemsBase.saturation that automatically forms the equivalent to nonlinearity(x->clamp(x, -0.7, 0.7)) while at the same time making sure the function has a recognizable name when the system is printed
using ControlSystemsBase: saturation
saturation(0.7)ControlSystemsBase.HammersteinWienerSystem{Float64}
P: StateSpace{Continuous, Float64}
D =
0.0 1.0
1.0 0.0
Continuous-time state-space model
Nonlinearities: Function[saturation(0.7)]See also ControlSystemsBase.ratelimit that saturates the derivative of a signal.
Non-zero operating point
It's common to linearize nonlinear systems around some operating point. We may make use of the helper constructor ControlSystemsBase.offset to create affine functions at the inputs and outputs of the linearized system to, e.g.,
- Make sure that simulations result are given in the original coordinates rather than in the coordinates of the linearization point.
- Allow nonlinearities that are added back after the linearization (such as saturations) to operate with their original parameters.
We will demonstrate a composite usage of offset and saturation below. The system we'll consider is a linearized model of a quadruple-tank process;
The system is linearized around the operating point
xr = [10, 10, 4.9, 4.9] # reference state
ur = [0.263, 0.263] # control input at the operating pointand is given by
using LinearAlgebra
kc, k1, k2, g = 0.5, 1.6, 1.6, 9.81
A1 = A3 = A2 = A4 = 4.9
a1, a3, a2, a4 = 0.03, 0.03, 0.03, 0.03
h01, h02, h03, h04 = xr
T1, T2 = (A1/a1)sqrt(2*h01/g), (A2/a2)sqrt(2*h02/g)
T3, T4 = (A3/a3)sqrt(2*h03/g), (A4/a4)sqrt(2*h04/g)
c1, c2 = (T1*k1*kc/A1), (T2*k2*kc/A2)
γ1, γ2 = 0.3, 0.3
# Define the process dynamics
A = [-1/T1 0 A3/(A1*T3) 0
0 -1/T2 0 A4/(A2*T4)
0 0 -1/T3 0
0 0 0 -1/T4]
B = [γ1*k1/A1 0
0 γ2*k2/A2
0 (1-γ2)k2/A3
(1-γ1)k1/A4 0 ]
C = kc*[I(2) 0*I(2)] # Measure the first two tank levels
D = 0
G = ss(A,B,C,D)A PID controller with a filter is given by
F = tf(1, [0.63, 1.12, 1])
Cpid = pid(0.26, 0.001, 15.9, form=:parallel)*F |> ssand to make the controller MIMO, we add a static pre-compensator that decouples the system at the the zero frequency.
iG0 = dcgain(G)
iG0 ./= maximum(abs, iG0)
C = (Cpid .* I(2)) * iG0The pumps (there are two of them) that service the tanks can only add liquid to the tanks, not remove liquid. The pump is thus saturated from below at 0, and from above at the maximum pump capacity 0.4.
using ControlSystemsBase: offset
umin = [0.0, 0.0]
umax = [0.4, 0.4]
yr = G.C*xr # Reference output
Gop = offset(yr) * G * offset(-ur) # Make the plant operate in Δ-coordinates
C_sat = saturation(umin, umax) * C # while the controller and the saturation operate in the original coordinatesControlSystemsBase.HammersteinWienerSystem{Float64}
P: StateSpace{Continuous, Float64}
A =
0.0 0.015625 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0
0.0 -1.5873015873015872 -1.777777777777778 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.015625 0.0
0.0 0.0 0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0 -1.5873015873015872 -1.777777777777778
B =
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
1.7142857142857144 4.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
4.0 1.7142857142857144 0.0 0.0
C =
0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0
0.025396825396825397 0.10317460317460318 6.309523809523809 0.0 0.0 0.0
0.0 0.0 0.0 0.025396825396825397 0.10317460317460318 6.309523809523809
D =
0.0 0.0 1.0 0.0
0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
Continuous-time state-space model
Nonlinearities: Function[saturation(0.0, 0.4), saturation(0.0, 0.4)]We now simulate the closed-loop system, the initial state of the plant is adjusted with the operating point x0-xr since the plant operates in Δ-coordinates
x0 = [2, 1, 8, 3] # Initial tank levels
plot(
plot(lsim(feedback(Gop*C_sat), yr, 0:1:3000, x0=[x0-xr; zeros(C.nx)]), layout=1, sp=1, title="Outputs", ylabel=""),
plot(lsim(feedback(C_sat, Gop), yr, 0:1:3000, x0=[zeros(C.nx); x0-xr]), layout=1, sp=1, title="Control signals", ylabel="")
)
hline!([yr[1]], label="Reference", l=:dash, sp=1, c=1)Duffing oscillator
In this example, we'll model and control the nonlinear system
\[\ddot x = -kx - k_3 x^3 - c \dot{x} + 10u\]
To do this, we first draw the block diagram
10u ┌───┐
──────►│+ │ ┌───┐ ┌───┐
┌────►│- │ ẍ │ 1 │ ẋ │ 1 │ x
│ ┌──►│- ├──►│ - ├┬─►│ - ├─┬──►
│ │ ┌►│- │ │ s ││ │ s │ │
│ │ │ └───┘ └───┘│ └───┘ │
│ │ │ │ │
│ │ │ ┌───┐ │ │
│ │ └───┤ c │◄─────┘ │
│ │ └───┘ │
│ │ │
│ │ ┌───┐ │
│ └─────┤ k │◄──────────────┤
│ └───┘ │
│ │
│ ┌───┐ ┌───┐ │
└───────┤ k³│◄──┤ x³│◄──────┘
└───┘ └───┘We see that the input $u$ passes through the inner velocity loop before reaching the output $x$, we can form this inner closed-loop transfer function using feedback(1/s, c), i.e., close the loop over an integrator by $-c$. This inner loop is then connected in series with another integrator an feedback loop is closed with $k_3 x^3 + kx =$ pos_loop_feedback in the feedback path. Notice how we multiply the final system with 10 from the right to get the input gain correct, for nonlinear systems, 10*sys and sys*10 are not always equivalent!
using ControlSystems, Plots
using ControlSystemsBase: nonlinearity
k = 10
k3 = 2
c = 1
s = tf("s")
cube = nonlinearity(x->x^3)
vel_loop = feedback(1/s, c)
pos_loop_feedback = (k3*cube + k)
duffing = feedback(vel_loop/s, pos_loop_feedback)*10
plot(step(duffing, 20), title="Duffing oscillator open-loop step response")We now show how we can make use of the circle criterion to prove stability of the closed loop. The function circle_criterion below plots the Nyquist curve of the loop-transfer function and figures out the circle to avoid by finding sector bounds for the static nonlinearity $f(x) = x^3$. We then choose a controller and check that it stays outside of the circle. To find the sector bounds, we choose a domain to evaluate the nonlinearity over. The function $f(x) = x^3$ goes to infinity faster than any linear function, and the upper sector bound is thus ∞, but if we restrict the nonlinearity to a smaller domain, we get a finite sector bound:
function circle_criterion(L::ControlSystemsBase.HammersteinWienerSystem, domain::Tuple; N=10000)
fun = x->L.f[](x)/x
x = range(domain[1], stop=domain[2], length=N)
0 ∈ x && (x = filter(!=(0), x)) # We cannot divide by zero
k1, k2 = extrema(fun, x)
f1 = plot(L.f[], domain[1], domain[2], title="Nonlinearity", lab="f(x)", xlab="x")
plot!(x, [k1.*x k2.*x], lab=["k1 = $(round(k1, sigdigits=2))" "k2 = $(round(k2, sigdigits=2))"], l=(:dash), legend=:bottomright)
p1 = -1/k2 # Close to origin
p2 = -1/k1 # Far from origin
c = (p1 + p2)/2
r = (p2 - p1)/2
Lnominal = sminreal(ss(L.A, L.B1, L.C1, L.D11, L.P.timeevol))
f2 = nyquistplot(Lnominal)
if p2 < -1000 # Due to bug in plots
vspan!([-1000, p1], fillalpha=0.7, c=:red, primary=false)
else
th = 0:0.01:2pi
Cs,Ss = cos.(th), sin.(th)
plot!(r.*Cs .+ c, r.*Ss, fill=true, fillalpha=0.7, c=:red, primary=false)
end
plot(f1,f2)
end
C = pid(2, 0, 1, form=:parallel)*tf(1, [0.01,1])
f1 = circle_criterion(duffing*C, (-1, 1))
plot!(sp=2, ylims=(-10, 3), xlims=(-5, 11))
f2 = plot(step(feedback(duffing, C), 8), plotx=true, plot_title="Controlled oscillator disturbance step response", layout=4)
plot(f1,f2, size=(1300,800))