# 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 C
```

```
ControlSystemsBase.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 point
```

and 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 |> ss
```

and 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)) * iG0
```

The 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 coordinates
```

```
ControlSystemsBase.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)
```

The state vector resulting from the call to `feedback`

is comprised of the concatenated states of the first and second arguments, i.e., `feedback(C_sat, Gop)`

has the state vector `[C_sat.x; Gop.x]`

while `feedback(Gop*C_sat)`

has the state vector of `Gop*C_sat`

which is starting with the first operand, `[Gop.x; C_sat.x]`

.

### 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))
```

Since we evaluated the nonlinearity over a small domain, we should convince ourselves that we indeed never risk leaving this domain.

In the example above, the circle turns into a half plane since the lower sector bound is 0. The example below chooses another nonlinearity

\[f(x) = x + \sin(x)\]

to get an actual circle in the Nyquist plane.

```
wiggly = nonlinearity(x->x+sin(x)) # This function is a bit wiggly
vel_loop = feedback(1/s, c)
pos_loop_feedback = (k3*wiggly + k)
duffing = feedback(vel_loop/s, pos_loop_feedback)*10
C = pid(2, 5, 1, form=:parallel)*tf(1,[0.1, 1])
f1 = circle_criterion(duffing*C, (-2pi, 2pi))
plot!(sp=2, ylims=(-5, 2), xlims=(-2.1, 0.1))
f2 = plot(step(feedback(duffing, C), 8), plotx=true, plot_title="Controlled wiggly oscillator disturbance step response", layout=5)
plot(f1,f2, size=(1300,800))
```

## Limitations

- Remember, this functionality is experimental and subject to breakage.
- Currently only
`Continuous`

systems supported. - No nonlinear root-finding is performed during simulation. This limits the kinds of systems that can be simulated somewhat, in particular, no algebraic loops are allowed.
- A lot of functions that expect linear systems will not work for nonlinear systems (naturally).

## Possible future work

- Discrete-time support.
- Basic support for nonlinear analysis such as stability proof through the circle criterion etc. In particular, predefined nonlinear functions may specify sector bounds for the gain, required by the circle-criterion calculations.
- Additional nonlinear components, such as
- Integrator anti-windup
- Friction models

## See also

More advanced nonlinear modeling is facilitated by ModelingToolkit.jl (MTK) and ModelingToolkitStandardLibrary.jl. The tutorials

- Modeling for control using ModelingToolkit
- Disturbance modeling in ModelingToolkit
- Modal analysis of a series of masses and springs using MTK

show how to use these packages to model and simulate control systems.

## Docstrings

`ControlSystemsBase.nonlinearity`

— Function```
nonlinearity(f)
nonlinearity(T, f)
```

Create a pure nonlinearity. `f`

is assumed to be a static (no memory) nonlinear function from $f : R -> R$.

The type `T`

defaults to `Float64`

.

NOTE: The nonlinear functionality in ControlSystemsBase.jl is currently experimental and subject to breaking changes not respecting semantic versioning. Use at your own risk.

**Example:**

Create a LTI system with a static input nonlinearity that saturates the input to [-1,1].

`tf(1, [1, 1])*nonlinearity(x->clamp(x, -1, 1))`

See also predefined nonlinearities `saturation`

, `offset`

.

Note: when composing linear systems with nonlinearities, it's often important to handle operating points correctly. See `ControlSystemsBase.offset`

for handling operating points.

`ControlSystemsBase.offset`

— Function`offset(val)`

Create a constant-offset nonlinearity `x -> x + val`

.

NOTE: The nonlinear functionality in ControlSystemsBase.jl is currently experimental and subject to breaking changes not respecting semantic versioning. Use at your own risk.

**Example:**

To create a linear system that operates around operating point `y₀, u₀`

, use

`offset_sys = offset(y₀) * sys * offset(-u₀)`

note the sign on the offset `u₀`

. This ensures that `sys`

operates in the coordinates `Δu = u-u₀, Δy = y-y₀`

and the inputs and outputs to the offset system are in their non-offset coordinate system. If the system is linearized around `x₀`

, `y₀`

is given by `C*x₀`

. Additional information and an example is available here https://juliacontrol.github.io/ControlSystemsBase.jl/latest/lib/nonlinear/#Non-zero-operating-point

`ControlSystemsBase.saturation`

— Function```
saturation(val)
saturation(lower, upper)
```

Create a saturating nonlinearity. Connect it to the output of a controller `C`

using

`Csat = saturation(val) * C`

```
y▲ ────── upper
│ /
│ /
│/
──────────┼────────► u
/│
/ │
/ │
lower────
```

NOTE: The nonlinear functionality in ControlSystemsBase.jl is currently experimental and subject to breaking changes not respecting semantic versioning. Use at your own risk.

Note: when composing linear systems with nonlinearities, it's often important to handle operating points correctly. See `ControlSystemsBase.offset`

for handling operating points.

`ControlSystemsBase.ratelimit`

— Function```
ratelimit(val; Tf)
ratelimit(lower, upper; Tf)
```

Create a nonlinearity that limits the rate of change of a signal, roughly equivalent to $1/s ∘ sat ∘ s$. `Tf`

controls the filter time constant on the derivative used to calculate the rate. NOTE: The nonlinear functionality in ControlSystemsBase.jl is currently experimental and subject to breaking changes not respecting semantic versioning. Use at your own risk.

`ControlSystemsBase.deadzone`

— Function```
deadzone(val)
deadzone(lower, upper)
```

Create a dead-zone nonlinearity.

```
y▲
│ /
│ /
lower │ /
─────|──┼──|───────► u
/ │ upper
/ │
/ │
```

Note: when composing linear systems with nonlinearities, it's often important to handle operating points correctly. See `ControlSystemsBase.offset`

for handling operating points.

`ControlSystemsBase.linearize`

— Function`linearize(sys::HammersteinWienerSystem, Δy)`

Linearize the nonlinear system `sys`

around the operating point implied by the specified Δy

```
┌─────────┐
y◄───┤ │◄────u
│ P │
Δy┌───┤ │◄───┐Δu
│ └─────────┘ │
│ │
│ ┌───┐ │
│ │ │ │
└─────►│ f ├───────┘
│ │
└───┘
```

`A, B = linearize(f, x, u, args...)`

Linearize dynamics $ẋ = f(x, u, args...)$ around operating point $(x,u,args...)$ using ForwardDiff. `args`

can be empty, or contain, e.g., parameters and time `(p, t)`

like in the SciML interface. This function can also be used to linearize an output equation `C, D = linearize(h, x, u, args...)`

.