In my group, many use Matlab and Simulink for control engineering. It is very easy for them since this his how you learn control engineering, by block diagrams. Regardless of the user interface, what options would I have to do this in Julia? Specifically simulation of a dynamical system (plant) and some form of causal controller (MPC, PID, custom rule-based logic, observers, sliding mode control, …).
There is:
ModelingToolkit.jl, for the dynamical system, but can it do the causal part? The discret system is still experimental and I am not sure if I can implement any function u = C(y,x,p,t). I cannot have it be symbolic and also it often needs some kind of state.
Causal.jl, which seems to be the closest to Simulink but sadly unmaintained (the author finished his PhD ). To me it seems perfect for loosely coupled co-simulation, since the right hand side can just be any function. One block could even be an MTK model (the plant).
DifferentialEquations.jl, which is the most low-level. Hybrid systems seem to be a bit cumbersome since the discrete part basically has to live in the callback of the dynamic part. (As mentioned in this 6 year old post DiffEqs: Hybrid (continuous+discrete) system / periodic callback)
My question now. Is there anything complete enough to bother my colleagues with?
As a use case, consider the classic text book problem of a two tank system and different control approaches for that. But of course it should handle larger systems.
Hybrid control is not yet very well supported by MTK. We just use MTK to create the model, and then use the integrator interface of DifferentialEquations to run the continues model timestep per timestep, the controller is a discrete time controller and calculates the new plant inputs on each timestep.
Please ask more specific questions to explain what you are missing.
Please, be aware that Dyad is a commercial software, free for academic use. You can do the same without Dyad, using MIT licensed packages only, but then you have less nice tutorials and a little bit extra work.
While creating a basic controller in Simulink is faster, doing the same with Julia is satisfying: I translated a Simulink model of a wind turbine + estimator + controller from Simulink to Julia, and the code was:
much more readable, mainly 35 differential algebraic equations
about 5000 times faster
the linearisations around different operating points much more accurate
A ready to use Simulink replacement that I can propose my colleagues. It is not about system analysis and controller synthesis (pole placement, transfer functions,…). It is about being able to craft hybrid systems and simulate them, without worrying about callbacks and stuff. My understanding is that MTK would be the Simscape, Modelica alternative. But what about the rest of the system (the causal part)? It might be feasible to make an FMI with MTK/Dyad and use it in Simulink.
Unfortunately, the expensive part in engineering is the engineer.
I just write it as Julia code. Example of a dual-loop controller:
@with_kw mutable struct ParkingControllerSettings @deftype Float64
dt
# turn rate controller settings
kp_tr=0.06
ki_tr=0.0012
# outer controller (heading/ course) settings
kp=15
ki=0.5
# NDI block settings
va_min = 5.0 # minimum apparent wind speed
va_max = 100.0 # maximum apparent wind speed
k_ds = 2.0 # influence of the depower settings on the steering sensitivity
c1 = 0.048 # v9 kite model
c2 = 0 # a value other than zero creates more problems than it solves
end
mutable struct ParkingController
pcs::ParkingControllerSettings
pid_tr::DiscretePID
pid_outer::DiscretePID
last_heading::Float64
chi_set::Float64
last_ndi_gain::Float64
end
function ParkingController(pcs::ParkingControllerSettings; last_heading = 0.0)
Ts = pcs.dt
pid_tr = DiscretePID(;K=pcs.kp_tr, Ti=pcs.kp_tr/ pcs.ki_tr, Ts)
pid_outer = DiscretePID(;K=pcs.kp, Ti=pcs.kp/ pcs.ki, Ts)
return ParkingController(pcs, pid_tr, pid_outer, last_heading, 0, 0)
end
"""
linearize(pcs::ParkingControllerSettings, psi_dot, psi, elevation, v_app; ud_prime=0)
Nonlinear, dynamic inversion block (NDI) according to Eq. 6.4 and Eq. 6.12.
Parameters:
- psi_dot: desired turn rate in radians per second
- psi: heading in radians
- elevation: elevation angle in radians
- v_app: apparent wind speed in m/s
- ud_prime: depower setting in the range of 0 to 1, 0 means fully powered, 1 means fully depowered
"""
function linearize(pc::ParkingController, psi_dot, psi, elevation, v_app; ud_prime=0)
pcs = pc.pcs
# Eq. 6.13: calculate va_hat
va_hat = clamp(v_app, pcs.va_min, pcs.va_max)
# Eq. 6.12: calculate the steering from the desired turn rate
u_s = (1.0 + pcs.k_ds * ud_prime) / (pcs.c1 * va_hat) * (psi_dot - pcs.c2 / va_hat * sin(psi) * cos(elevation))
if abs(psi_dot) < 1e-6
ndi_gain = pc.last_ndi_gain
else
ndi_gain = clamp(u_s / psi_dot, -20, 20.0)
end
pc.last_ndi_gain = ndi_gain
return u_s, ndi_gain
end
"""
navigate(fpc, limit=50.0)
Calculate the desired flight direction chi_set using great circle navigation.
Limit delta_beta to the value of the parameter limit (in degrees).
"""
function navigate(pc::ParkingController, azimuth, elevation; limit=50.0)
phi_set = 0.0 # azimuth
beta_set = deg2rad(77) # zenith
beta = elevation
phi = azimuth
# println("beta: $(rad2deg(beta)), phi: $(rad2deg(phi))")
r_limit = deg2rad(limit)
if beta_set - beta > r_limit
beta_set = beta + r_limit
elseif beta_set - beta < -r_limit
beta_set = beta - r_limit
end
y = sin(phi_set - phi) * cos(beta_set)
x = cos(beta) * sin(beta_set) - sin(beta) * cos(beta_set) * cos(phi_set - phi)
pc.chi_set = atan(y, x)
end
"""
calc_steering(pc::ParkingController, heading; elevation=0.0, v_app=10.0, ud_prime=0.0)
Calculate rel_steering and ndi_gain from the actual heading, elevation, and apparent wind speed.
Parameters:
- pc: parking controller
- heading: actual heading in radians
- elevation: elevation angle in radians
- v_app: apparent wind speed in m/s
- ud_prime: depower setting in the range of 0 to 1, 0 means fully powered, 1 means fully depowered
"""
function calc_steering(pc::ParkingController, heading, chi_set; elevation=0.0, v_app=10.0, ud_prime=0.0)
# calculate the desired turn rate
heading = wrap2pi(heading) # a different wrap2pi function is needed that avoids any jumps
psi_dot_set = pc.pid_outer(wrap2pi(chi_set), heading)
psi_dot = (wrap2pi(heading - pc.last_heading)) / pc.pcs.dt
pc.last_heading = heading
psi_dot_in = pc.pid_tr(psi_dot_set, psi_dot)
# linearize the NDI block
u_s, ndi_gain = linearize(pc, psi_dot_in, heading, elevation, v_app; ud_prime)
u_s, ndi_gain, psi_dot, psi_dot_set
end
I find this pretty straightforward. But if you say that your colleagues cannot program, that they are used to GUI tools only, then there is no easy way to replace Simulink with Julia.
The closest you can get is Dyad, but it is very new and I am not sure if the GUI part is already released. And it is expensive (for commercial purposes).
Implementing the model with MTK and exporting it as FMI is also a very nice idea, but not very mature yet. Until now you can only use one instance of a Julia FMI model in Simulink.
By the way, you do not need any callbacks if you implement a controller as shown above and the model using MTK as long as you use the integrator interface.
I think this is not what he is saying. Anyway, whether or not the control engineer is a good programmer, the graphical language of block diagrams (in the Simulink style, that is, essentially just signal flow graphs) is just so much convenient for both writing and reading the rules for manipulating the signals, aka the control law. After all, that is why they have been used in industry.
The hybrid systems are not quite released yet but should beta later this year.
I’d maybe wait for that if the discrete part is needed. But if you want earlier access, we’d be happy to work with you and start getting examples and test cases.
What Dyad really needs to do to gain adoption is to get into universities first! As someone who recently graduated from electrical engineering, control is still all done in matlab. The courses use matlab, the professors practically only know matlab, and in the end projects and internships are also done in matlab. Fortunately, some professors are forward-looking and use Julia (that’s how I got introduced to it), but they are rare.
I follow MTK for some years now, and even contributed some small things, years ago. For me the issue is not the lack of a GUI, although it is nice to have, it is the state of the interface and a solid foundation.
Instantiating components, connecting them and simulating them is fine in code as well. It might even be better since you can programmatically change stuff. However, doing low-level, manual integrator stepping is not feasible. It contradicts the composabiliy argument. Isn’t that what Causal.jl already somewhat solved (without any symbolic simplification and stuff)?
For a broader adoption, documentation, examples and tutorials are the most important thing I think. I would regard MATLAB’s documentation as pretty substantial. You will find almost everything in many domains, with examples, images, API documentation and reference literature. And everything is publicly accessible, so that a potential customer can have his/her “this looks cool moment”. But of course it is unfair, MATLAB has gotten millions of man-hours and a lot of industry money into it.
For instance what I am missing in MTK, is documentation on systems that use stream connectors such as thermal fluid flow. They are implemented since 2021! but undocumented without any examples for some reason.
For me, a real world use case would be a district heating network with 200 customers, each with a digital controller and controllers for the heating plant. I know this is hard, and others have tried and failed, but this is Julia after all.
Here is a five year old example I made, demonstrating the wind-up effect. As you can see the syntax is ok, however, the right hand side of the PI-controller with back calculation looks awful since you need to get rid of algebraic loops.
And Causal.jl is more like loosely coupled co-simulation than what MTK does. But it might be enough for you.
using Causal, Plots
"""
System.
"""
@def_ode_system mutable struct Sys{RH,RO,ST,IP,OP} <: AbstractODESystem
righthandside::RH = function ode(dx, x, u, t)
dx[1] = x[2]
dx[2] = -1. * x[1] - 0.5 * x[2] + u[1](t)
end
readout::RO = (x, u, t) -> 0.9 * x[1] + x[2]
state::ST = [0., 0.]
input::IP = Inport()
output::OP = Outport()
end
"""
Time continuos PI controller.
"""
@def_ode_system mutable struct PIController{RH,RO,ST,IP,OP} <: AbstractODESystem
kp::Float64 = 1.
Tn::Float64 = 0.1
righthandside::RH = (dx, x, u, t, kp = kp, Tn = Tn) -> (dx[1] = kp / Tn * u[1](t))
readout::RO = (x, u, t, kp = kp) -> kp * u[1](t) + x[1]
state::ST = [0.]
input::IP = Inport()
output::OP = Outport()
end
"""
Time continuos PI controller with actuator saturation.
"""
@def_ode_system mutable struct PIControllerWithSat{RH,RO,ST,IP,OP} <: AbstractODESystem
kp::Float64 = 1.
Tn::Float64 = 0.1
lb::Float64 = -1.
ub::Float64 = 1.
righthandside::RH = (dx, x, u, t, kp = kp, Tn = Tn, lb = lb, ub = ub) -> (dx[1] = kp / Tn * u[1](t))
readout::RO = (x, u, t, kp = kp, ub = ub, lb = lb) -> max(min(kp * u[1](t) + x[1], ub), lb)
state::ST = [0.]
input::IP = Inport()
output::OP = Outport()
end
"""
Time continuos PI controller actuator saturation and anti-windup.
"""
@def_ode_system mutable struct PIControllerWithSatAndAntiWindup{RH,RO,ST,IP,OP} <: AbstractODESystem
kp::Float64 = 1.
Tn::Float64 = 0.1
lb::Float64 = -1.
ub::Float64 = 1.
Ta::Float64 = 0.1
righthandside::RH = function ode(dx, x, u, t, kp=kp, Tn=Tn, Ta=Ta, lb=lb, ub=ub)
dx[1] = kp / Tn * u[1](t) + 1 / Ta * (-(x[1] + kp * u[1](t)) + max(min(kp * u[1](t) + x[1], ub), lb))
end
readout::RO = (x, u, t, kp = kp, lb = lb, ub = ub) -> max(min(kp * u[1](t) + x[1], ub), lb)
state::ST = [0.]
input::IP = Inport()
output::OP = Outport()
end
# Construct the model
@defmodel model1 begin
@nodes begin
step = StepGenerator(delay=1.)
sys = Sys()
sum = Adder(signs=(+, -))
pi_contr = PIController(kp=3., Tn=0.5)
writer = Writer(input=Inport(3))
mem = Memory(delay=0.01)
end
@branches begin
step[1] => sum[1]
sys[1] => mem[1]
mem[1] => sum[2]
sum[1] => pi_contr[1]
pi_contr[1] => sys[1]
step[1] => writer[1]
sys[1] => writer[2]
pi_contr[1] => writer[3]
end
end
# Construct the model without anti-windup measure
@defmodel model2 begin
@nodes begin
step = StepGenerator(delay=1.)
sys = Sys()
sum = Adder(signs=(+, -))
pi_contr = PIControllerWithSat(kp=3., Tn=0.5, lb=-1.5, ub=1.5)
writer = Writer(input=Inport(3))
mem = Memory(delay=0.01)
end
@branches begin
step[1] => sum[1]
sys[1] => mem[1]
mem[1] => sum[2]
sum[1] => pi_contr[1]
pi_contr[1] => sys[1]
step[1] => writer[1]
sys[1] => writer[2]
pi_contr[1] => writer[3]
end
end
# Construct the model with anti-windup measure
@defmodel model3 begin
@nodes begin
step = StepGenerator(delay=1.)
sys = Sys()
sum = Adder(signs=(+, -))
pi_contr = PIControllerWithSatAndAntiWindup(kp=3., Tn=0.5, lb=-1.5, ub=1.5, Ta=0.1)
writer = Writer(input=Inport(3))
mem = Memory(delay=0.01)
end
@branches begin
step[1] => sum[1]
sys[1] => mem[1]
mem[1] => sum[2]
sum[1] => pi_contr[1]
pi_contr[1] => sys[1]
step[1] => writer[1]
sys[1] => writer[2]
pi_contr[1] => writer[3]
end
end
# call simulate
sim1 = simulate!(model1, 0, 0.01, 6.)
sim2 = simulate!(model2, 0, 0.01, 6.)
sim3 = simulate!(model3, 0, 0.01, 6.)
# Read and plot data
t, x = read(getnode(model1, :writer).component)
p1 = plot(t, x[:, 1], label="r(t)", xlabel="t", legend=:bottomright)
plot!(p1, t, x[:, 2], label="y(t) without actuator saturation", xlabel="t")
p2 = plot(t, x[:, 3], label="u(t) without actuator saturation", xlabel="t")
t, x = read(getnode(model2, :writer).component)
plot!(p1, t, x[:, 2], label="y(t) without anti-windup", xlabel="t")
plot!(p2, t, x[:, 3], label="u(t) without anti-windup", xlabel="t")
t, x = read(getnode(model3, :writer).component)
plot!(p1, t, x[:, 2], label="y(t) with anti-windup", xlabel="t")
plot!(p2, t, x[:, 3], label="u(t) with anti-windup", xlabel="t")
plot(p1, p2, layout=(2, 1))
MATLAB is a paid product, you pay for what you get. MTK is mostly built on community contributions, so we, the community, also have to put in the work. A lot of the basic stuff is still missing, for example a short while back I contributed the diode component. Based on the diode component I can now build single diode solar cells, which then can be used to build more comprehensive simulations. Comparing that to MATLAB, we are still far behind. But at the same time, I believe that the foundation is in place thanks to the work put in by people like Chris. But we, the users, also have to put in some work to make it happen.
). To me it seems perfect for loosely coupled co-simulation, since the right hand side can just be any function. One block could even be an MTK model (the plant).