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.
As a university lecturer (senior lecturer) in control systems theory over about 30 years, I was frustrated that I was expected to teach MATLAB/Simulink (which I did in laboratory classes) but that students couldn’t work on it at home and that the number of toolboxes we had was limited. I switched to Scilab in my notes (skipping the ever-so-slow Octave) so that students could experiment with it without downloading hacked and illegal versions of MATLAB. However, I felt sure others had the same problem, which is how in 2015/16 I stumbled on Julia and became an advocate.
I also found that students who depended on MATLAB/Simulink had not really learnt control theory except as a set of isolated and disjointed recipes. They couldn’t implement them in code because if you want to do actual programmatical discrete-time control, you’re faced with the problem MATLAB gives you in that you have no control over how it writes mex files. In industry, how do you support the code? Handling telemetry message queues at high speed (hello ZMQ.jl), multi-threading, co-routines etc, storing data in a fast time-series database QuestDB seems to be only possible with code and I found Julia perfect for the job. If you must interface with MATLAB (co-simulation), you can use FMI.jl.
I would advocate using block diagrams for modular design but do the implementation in code that you can optimise yourselves. I couldn’t see how to use MTK for realtime stuff with multiple streams of measurements and I felt too far removed from the problem and so preferred my own code. At least then I knew the errors where mine.
Good discussion. My thoughts always focus on solving the three language problem and not just the two language problem to get to the embedded controller.
I think key about the controller design is how is that embeddable into a real time controller. Being able to prototype the discrete portion in Julia is great, but ultimately if there isn’t a clear and verifiable path to running the same code on an embedded controller, it is just academic.
If you must do the two language prototype problem for an embedded controller with Julia being transcribed to another language, the simulation needs to be flexible enough to allow for verification and validation inside the simulation. Ideally Julia would autogenerate the c equivalent and verify it.
When this requires some of the logic being implemented in fixed point logic or VHDL type architectures, we start to see the significant gap between MATLAB and Julia which needs to be closed.
In the aviation industry only subsets of C, C++ and ADA are supported, but there are large areas of control system applications that are less restrictive.
MATLAB can generate code for such systems, but it costs a lot of $$$. Is that reasonable to expect from open source software? I think it would require an entire engineering department to build and maintain such tools.
I think there is an engineering team working on Dyad Studio, and as far as I know Boeing is one of their customers. So I guess this is being worked on. How much of it will be published with an MIT license we will see.
We have chosen to publish Dyad Studio under a Source Available license. This allows us to make the code publicly available on GitHub and people can use it for free for educational and non-commercial personal purposes, with licenses available for commercial use. In addition, Dyad’s standard libraries are BSD3 open source licensed to enable open bug reporting and rapid improvements by community contributions.
That’s great of course. But then I also expect Dyad to follow the same approach as outlined in the paper by JuliaHub, although Dyad is not mentioned explicitly. So no C codegen.
MTK performs code generation from the symbolic representation to Julia code (C-code generation for this step is also available).
I think skipping C code generation is the preferred option, but the C code generation for digital controllers might also come if there is customer demand for it.
). 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).