I have a discrete MIMO system in state space form and I’m trying to convert the internal feedback path from z^-1 to a continuous delay (i.e. exp(-s*T)).
In MATLAB, I have the fairly clunky but working solution to plug the system matrices into the right places of setDelayModel, but I can’t find any equivalent in Julia.
I need this because I need to analyse a discrete controller + continuous system from the continuous side without the phase lag caused by discretization. Only a delay can give me an exact match up to the nyquist frequency, ZOH transform, Padé etc. don’t do that.
I guess extracting the system matrices and connecting them as gains rather than individual systems should work, but I don’t even know how to incorporate a delay into that.
Two associated questions I could not find the answer to:
Do MIMO delay systems even exist in Julia beyond arrays of SISO systems?
Can I name their inputs and outputs?
All I need is to compute the frequency response, no need for stability analysis or time domain simulations.
I am not sure I understand exactly what you are trying to do here. In particular, I’m not sure what approach you are suggesting to convert the discrete-time controller to continuous time to “analyze from the continuous side”?
If you only want frequency responses, you can always perform computations with the frequency responses directly, i.e., compute the frequency response of the discrete-time part G_d(i\omega) and the continuous-time part G_c(i\omega) independently and perform frequency-wise computations with G_d(i\omega) and G_c(i\omega).
In ControlSystems.jl, you create a delay term using delay(L) and this can then be used like any ordinary LTI system. All system types in ControlSystems are inherently MIMO, including the DelayLTISystem type.
Thanks for your swift reply.
I have a continuous time plant with a discrete time controller. But the plant has an input (called ug) that is not directly manipulated by the controller. I need to assess how the discrete controller affects the transfer function between the continuous input ug and one of the (continuous) outputs.
The standard method of combining a discrete controller with a continuous plant is to just discretize the plant using ZOH and slap the discrete controller onto it. The problem with that method is that this will give me ZOH behavior and thus distorted frequency response for all inputs and outputs - even those that I’d like to know continuous time behavior for.
My solution to the problem is to use a continuous ZOH for only the plant input that the controller manipulates and model everything in the continuous domain. Because I need good (maximum) accuracy near the nyquist frequency, I can’t just approximate my discrete controller with d2c or similar, I need an exact match. That exact match is provided by replacing the delay (z^-1) by a continuous delay, i.e. exp(-s*T). The only difference being that the resulting transfer function has a life beyond the nyquist frequency that doesn’t represent reality, but I only care what the controller does to the system, so I only care about frequencies up to nyquist.
That delay systems don’t have names is sad, but I suppose I’m going to have to deal with that.
I see what problem you are trying to solve now. I’m still not sure I follow what you mean with
The continuous-time transfer function of the ZoH operator is given by
\dfrac{1 - e^{-sT}}{sT}
or in ControlSystems.jl language
s = tf('s')
Ts = 1 # Sample interval
Z = (1 - delay(Ts))/s
Just replacing the discrete-time controller with a transfer function that is a rational function of the delay element will give a system with continuously varying output?
The alternatives I can think of are
Perform the computations directly in the frequency domain like outlined above.
Use Tustin approximation with frequency prewarping to get the correct frequency response for higher frequencies.
Manually construct a transfer function in the delay(Ts)= e^{-sT_s} element, e.g.,
d = delay(Ts)
d^2 + 2d + 1
would model
e^{-2sT_s} + 2e^{-sT_s} + 1
Internally yes, but from the users perspective we do not expose any SISO system, we expose a single transfer-function type and this is MIMO, but is represented internally as a matrix of SISO rational functions. Do you have any particular reason to use transfer functions? They tend to be numerically worse in most ways when compared to statespace representations. Our transfer-function type does not support delays either, all delay systems in ControlSystems.jl are represented on LFT form where an array of pure delay elements is in feedback with a partitioned statespace system.
z^{-1} is a delay of the sampling interval T and has no rational correspondence in the s-domain, but an exact irrational equivalent: e^{-sT}. Sure, the output of the latter is continuous while the output of the former isn’t, but their frequency responses are exactly equal until the nyquist frequency. Using that in combination with the ZOH operator on one of my system inputs allows me to compute a frequency response without ZOH operator for the other system input, while still having an exact representation of my discretely controlled system (in the frequency domain - the time domain will be off because there is no discretization).
Building the whole system from transfer functions would be possible, but it’s unwieldy. I already have the named state space representations of my plant and my controller, having to convert and connecting between various unnamed inputs is complicated and error prone. But since I’m probably going to have to try that, how would I conveniently replace z^{-1} \rightarrow e^{-sT} in a given discrete time transfer function?
Tustin doesn’t work, I need accuracy across the entire frequency range - I also wouldn’t be satisfied with an approximation because an exact computation is absolutely feasible.
I don’t think there is a very conventient way of doing this, it’s not really a use case we have considered yet.
Something like this perhaps?
function d2dc(sys::TransferFunction{<:Discrete})
ControlSystemsBase.issiso(sys) || throw(ArgumentError("Only SISO systems are supported"))
T = sys.Ts
n, d = reverse(numvec(sys)[]), reverse(denvec(sys)[])
ne = [n[i]*(i == 1 ? 1 : delay(-T*(i-1))) for i in 1:length(n)] |> sum
de = [d[i]*(i == 1 ? 1 : delay(-T*(i-1))) for i in 1:length(d)] |> sum
ne * feedback(1, de-1)
end
Ts = 1 # Sample interval
P = tf(0.1, [1, 0.1, 0.1])
Pd = c2d(P, Ts) # A discrete-time system
w = exp10.(LinRange(-2, 1, 200))
Pdc = d2dc(Pd)
bodeplot(Pd, w, lab="Discrete"); bodeplot!(Pdc, w, lab="Continuous with delays")
I have temporarily given up and reverted to MATLAB, though I am determined to pick up Julia again very soon.
A little background for what I’m doing:
I’m trying to evaluate a power electronic converter (“inverter”) control in terms of passivity. A power system component is “passive” at a particular frequency if the real part of its impedance (or admittance, the transfer function from grid voltage to grid current) is positive at that particular frequency. That is, it’s resistance is positive - a negative resistance is a potential source of instability. If every component of a power system is passive at a given frequency, the entire system will be stable at that frequency. This is useful, because I no longer have to analyze the entire system, which often isn’t feasible. If my component is passive across the entire frequency range, I can guarantee it will not cause problems when connected to the larger power system.
The problem with this analysis is that the converter hardware is continuous, interaction with the grid (this is where I need to calculate the impedance) is continuous, but the controller is discrete. A discrete controller can make a major difference for higher frequencies, so approximating it as continuous won’t give me accurate results at higher frequencies.
While a lot of people won’t bother with such details, I don’t think there’s a particular reason not to, because the computational effort shouldn’t be significantly higher than that of inaccurate approximations. Switching frequencies of around 3 kHz aren’t rare and neither are resonances at those frequencies - enough reason to make the effort and not have any bad surprises.
Here’s the corresponding one for statespace systems
using LinearAlgebra
function d2dc(sys::AbstractStateSpace{<:Discrete})
T = sys.Ts
A,B,C,D = ssdata(sys)
z = delay(-T)
LR = append([z for _ in 1:sys.nx]...) - ss(A + I)
C * feedback(I(sys.nx), LR)*B + D
end
sys = ssrand(2, 3, 2, Ts = 1) # A random discrete-time MIMO system
sysc = d2dc(sys)
bodeplot(sys, w, lab="Discrete SS")
bodeplot!(sysc, w, lab="Continuous SS with delays", l=:dash, size=(800, 800))
Cool, thanks for the detail! I was considering asking what your application was that required such precision in the frequency response close to the Nyquist frequency.
Wow, that’s amazing! Thank you for that quick implementation.
As for the tutorial page, it pretty much sums things up very well, I have a few points of criticism:
You wrote d2d_exact once instead of d2c_exact
It should be noted that a discrete controller doesn’t really have a frequency response above the nyquist frequency (it would cause aliasing for any frequencies above nyquist). So even with d2c_exact, care must be taken when interpreting results for those frequencies. For my use case, I know the system is passive above nyquist, so I just need results up to nyquist. I also noticed that bodeplots for discrete system extend beyond nyquist, which is questionable - that’s another topic though.
“If we choose option 1. in this case, we get a discretized…” → “If we choose option 1, we get a discretized…”
Is the function name temporary? I suppose it would be more intuitive if it became part of d2c, with a :delay option for example. I also think “delay” or the like is a better term than “exact”, because it describes what it does rather than what it is trying to achieve - the latter may not be true every time (see point 2 above).
I’m also not too sure how to model mixed sampling frequency systems. I suppose that’s for now beyond the scope of the tutorial page. I’m modelling the controller component with (much) higher sampling frequency as part of the continuous system and only the slower sampling frequency part as discrete, but that probably doesn’t capture the full behavior.
I have by the way managed to implement my model in Julia using the second method you wrote so I can perform analysis with it now. Doing that without being able to name inputs and outputs was a bit of a feat, so I’m looking forward to a more general named_system method
I will try out your implementation of d2c_exact tomorrow!
Thanks a lot for taking a look and for the feedback
Yeah, a Bode plot is perhaps not the best way of visualizing what happens above the Nyquist frequency. I believe that the magnitude computed above the Nyquist frequency will be correct for the fundamental component if we interpret the discrete system as having a sample operator at the input, at least if I read Ch 7.7 of Åström, Computer Controlled Systems, correctly. The frequency modulation introduced by sampling would still have to be considered separately though. I will insert a reference to this chapter in the tutorial and emphasize the caveats to the interpretation of the frequency response. I also forgot to include the transfer function of the hold operator in some of my examples.
It would be nice if it could be part of d2c as you say, but it isn’t for Julia-technical reasons. If passed a StateSpace model, d2c always returns another StateSpace model. d2c_exact returns a different type DelayLTISystem, and incorporating this feature into the same function would make d2ctype unstable with poor performance and compile-time regressions as result. The compiler would in this case not be able to infer the type of the output based on the types of the inputs to d2c.
d2c_exact works for me. Thanks again! I’m just wondering, are named non-state-space-systems far ahead on the roadmap? I used matrices of ones and zeros around my systems to match input to outputs, but that’s rather tedious.
I don’t know how Julia represents DelayLTISystems internally. In MATLAB it seems to be a collection of matrices that connect between integrators and delays. So a state space systems are just a subset of delay systems (see setDelayModel).
We represent them as an LFT between a partitioned statespace system and a diagonal operator of pure delays.
There is not really a roadmap ControlSystems.jl is an open-source project and has largely been developed through volunteer efforts. Nowadays, I maintain it partly on my free time with the help of occasional contributors, and partly during my work on JuliaSim. I estimate that implementing a more generic named type is going to be a significant effort, and even though I have plans for how such a type could look, incorporating not only delays, but things like static nonlinearities, uncertainties, time-varying components etc. as well, I am not sure when opportunity to spend the effort implementing this will arise.
I see, looks like delay DelayLTISystem looks a lot like what MATLAB also uses.
a generic named type would be awesome, it should enable constructing more complex systems without losing track of inputs and outputs. But of course, I understand that things may take time
