Feedback and inverses of MIMO transfer functions

[Krebbekx]

Hey all! I am not sure if this is the right place to request a feature. I am working on a Julia project that uses MIMO transfer functions as input. It is often handy to be able to construct these MIMO transfer functions by doing things like P = feedback(G_1, G_2) or even more simply G = P / Q.

When I try it with the current state of affairs, I get the error message MIMO TransferFunction feedback isn't implemented. or MIMO TransferFunction inversion isn't implemented yet. Is there any ambition to implement this soon?

[baggepinnen]

Hi there 👋

Computations involving transfer functions are numerically poor, so I usually recommend people to convert their transfer functions to state-space representation as soon as possible for better performance. Is there a reason not to do this in your case? ss(G_1) converts G_1 to StateSpace using a canonical basis on which numerical balancing is performed by default, making it behave well in computations.

[baggepinnen]

This PR does the conversion to ss automatically, performs the computation and then converts back to tf. You'll also get a warning indicating that it's preferable to convert to ss manually before the computation in order to avoid the conversion back to tf of the result.

• handle MIMO tf feedback and division by converting to ss #1007

[Krebbekx]

Thanks for the fast reply! I use the matrices per frequency in my algorithm, so for that I do really need the tf and not the ss representation of the system.

Another reason to raise the issue is that I try to promote Julia to my colleagues, who are mostly MATLAB users. Therefore, these basic operations such as feedback should not yield an error that gives the impression that the toolbox is very immature :).

I think your proposed solution is very good, since it concludes the explanation that ss is a better choice for numerical performance.

[Krebbekx]

Sure! I am using the TF of the plant G to find a controller K for the system. There may be various performance constraints on the sensitivity (1+GK)^{-1} and open-loop GK, and an optimization algorithm searches for the controller that satisfies the constraints.

The primal reason to use TF instead of SS representation is that, when working with LTI systems, you can swap the TF with frequency-response measurement data directly. Of course you can obtain an SS model via a system identification step, but the precise aim of the project is to skip the identification step.

[baggepinnen]

I don't understand what a TF provides you with here that a statespace system doesn't also provide. A statespace system implies a transfer function

$$C(sI - A)^{-1}B + D = \dfrac{N(s)}{D(s)}$$

and there is no problem computing a frequency response from a statespace system.

I perform similar optimization of controllers with frequency-domain constraints on sensitivity functions here
https://juliacontrol.github.io/ControlSystems.jl/stable/examples/automatic_differentiation/#Optimization-based-tuning%E2%80%93PID-controller
note how all system types involved in these examples are on statespace form.

[Krebbekx]

You are correct in stating that if you have an SS representation, it implies a TF representation.

In our case, we have a setup in the lab, and we perform FRF data measurements. These measurements, which are frequency-gain/phase tuples, are then used as transfer functions that are sampled at the said frequency points. I am therefore only interested in evaluating the transfer function of the controller at those sampled frequency points and multiplying that matrix of complex values with the FRF measurement data.

To remedy large condition numbers in the MIMO matrices, we use i/o normalization which is standard in H-infinity synthesis.

To use the SS framework, I have to identify an SS model from the FRF data, and the point of the research project is to skip this step. Of course, we will compare the approaches.

I also see that you use direct simulation in your cost function. This is something we do not do, which would indeed result in unreliable results when using TFs. We only optimize the controller based on appropriate gain bounds.

[baggepinnen]

I'm sorry, but I still don't see any issue with using a state-space representation for the controller, you can evaluate this at any frequency of your liking:

julia> Kss = ssrand(2,2,2)
StateSpace{Continuous, Float64}
A = 
 -1.8246362584299034  -0.4464801771538347
 -0.8662699700715216  -0.47685349079010264
B = 
 -0.9817875283388285  0.6455110870093906
  0.6397479473058993  0.6112225432150674
C = 
  1.7209912299136276  -0.3496797962362475
 -0.7519685741931539  -1.4022134124956138
D = 
 0.05506433091567787   0.48594341052564577
 0.44028274971774095  -0.07127874333997815

Continuous-time state-space model

julia> Ktf = tf(Kss)
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
Input 1 to output 1
0.05506433091567787s^2 - 1.7866246646160948s - 1.976260137879066
----------------------------------------------------------------
       1.0s^2 + 2.301489749220006s + 0.48331179955391124

Input 1 to output 2
0.44028274971774095s^2 + 0.8545164507504781s - 2.049756263119376
----------------------------------------------------------------
       1.0s^2 + 2.301489749220006s + 0.48331179955391124

Input 2 to output 1
0.48594341052564577s^2 + 2.0155805232145108s + 0.10050427049431673
------------------------------------------------------------------
        1.0s^2 + 2.301489749220006s + 0.48331179955391124

Input 2 to output 2
-0.07127874333997815s^2 - 1.5065157969744154s - 0.840436729358966
-----------------------------------------------------------------
        1.0s^2 + 2.301489749220006s + 0.48331179955391124

Continuous-time transfer function model

julia> special_frequency = 3.14
3.14

julia> freqresp(Ktf, special_frequency)
2×2 Matrix{ComplexF64}:
 -0.120745+0.505255im   0.640208-0.181559im
  0.565951+0.150034im  -0.234729+0.323598im

julia> freqresp(Kss, special_frequency)
2×2 Matrix{ComplexF64}:
 -0.120745+0.505255im   0.640208-0.181559im
  0.565951+0.150034im  -0.234729+0.323598im

Here's a complete example that mixes a statespace representation of a controller with a FrequencyResponseData object from ControlSystemIdentification.jl

using ControlSystemsBase, ControlSystemIdentification
G = tf(1, [1, 2])
w = exp10.(LinRange(-2, 2, 20))
Gf = FRD(w, G) # This creates a frequency-response data object

K = ssrand(1,1,2);

You can perform operations with this statespace controller and the FRD object, which automatically samples K at the same frequencies as those present in the FRD:

julia> K*Gf
Frequency (rad/s)
----------------
20-element Vector{Float64}:
   0.01
   0.01623776739188721
   0.026366508987303583
   0.04281332398719396
...
  14.38449888287663
  23.357214690901213
  37.926901907322495
  61.584821106602604
 100.0

Response
--------
20-element Vector{ComplexF64}:
    -2.7246366849192936 + 0.4074806700287081im
     -2.643252295593545 + 0.64503450928045im
    -2.4468634404585807 + 0.9822511002106853im
     -2.028689892707133 + 1.3696726253048837im
...
   0.005248368734117695 - 0.031358078055101705im
  0.0020105162106898006 - 0.019559403436560964im
  0.0007654451356195358 - 0.012104328981243383im
 0.00029073199220699866 - 0.007468220274811326im
 0.00011032653215554974 - 0.004602518317135966im

julia> feedback(K, Gf)
Frequency (rad/s)
----------------
20-element Vector{Float64}:
   0.01
   0.01623776739188721
   0.026366508987303583
   0.04281332398719396
...
  14.38449888287663
  23.357214690901213
  37.926901907322495
  61.584821106602604
 100.0

Response
--------
20-element Vector{ComplexF64}:
  3.0970529832214067 + 0.2749995147711004im
   3.051240447131617 + 0.4387682182575588im
   2.937740906980385 + 0.6812061573518066im
  2.6811858302947296 + 0.9914130092815417im
...
 0.46013269107459065 + 0.01038395373251116im
 0.46032881321268965 + 0.006404212246941321im
  0.4604035276098903 + 0.003946196880164585im

[Krebbekx]

Okay now I see your point. Indeed, the controller can perfectly fine be in an SS representation, since the only thing required are frequency domain samples. I just never want to compute a SS representation of the plant at any point.