# Trouble with ControlSystems.jl

**URL:** <https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783>\
**Category:** Modelling & Simulations\
**Tags:** control\
**Created:** [February 2, 2018, 5:37pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783 "2018-02-02T17:37:18Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![bo-lo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bo-lo/32/3317_2.png) [@bo-lo](https://discourse.julialang.org/u/bo-lo)\
**Post date:** [February 2, 2018, 5:37pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/1 "2018-02-02T17:37:18Z")

</div>

Hi

I am trying to use ControlSystems.jl to model circuits involving inductors, capacitors and resistances in the s domain.

## On my first try, I am getting a Zero denominator issue that I cannot debug. Here is the code (Julia 0.5.2):

```julia-auto
Pkg.update()
Pkg.add("ControlSystems")
using ControlSystems

rl = 0.2

L = 150.0e-9
dcr = 2.0e-3 # actual dcr is 0.2 mohm. I am adding the scaled top/bottom fet resistances.

nbulk = 12
cbulk = 470.0e-6
esr_bulk = 4.5e-3

ncer = 96
ccer = 14.0e-6 #22uF derated for 1V DC.
esr_cer = 2.0e-3

nph = 1

s = tf("s")

l_imp = dcr + L*s
bulk_imp = (esr_bulk/nbulk) + 1.0/(cbulk*nbulk*s)
cer_imp = (esr_cer/ncer) + 1.0/(ccer*ncer*s)
cer_load_imp = cer_imp*rl/(cer_imp + rl) #ceramic + load resistance.
den = (bulk_imp + cer_load_imp)
num = (bulk_imp * cer_load_imp)
tot_load_imp = num/den

```

* * *

Any help would be appreciated.

Thanks,  
Venkat.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 2, 2018, 11:39pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/2 "2018-02-02T23:39:30Z")

</div>

Please have a look at: [PSA: how to quote code with backticks](https://discourse.julialang.org/t/psa-how-to-quote-code-with-backticks/7530)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 3, 2018, 12:21am UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/3 "2018-02-03T00:21:26Z")

</div>

The code you pasted does not execute so it’s unfortunately not very easy to know how to help you.

---

<div class="post-metadata">

**Author:** ![bo-lo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bo-lo/32/3317_2.png) [@bo-lo](https://discourse.julialang.org/u/bo-lo)\
**Post date:** [February 3, 2018, 12:31am UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/4 "2018-02-03T00:31:56Z")

</div>

Thanks,  
I have updated my post.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [February 3, 2018, 5:07am UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/5 "2018-02-03T05:07:23Z")

</div>

On 0.6.2, but should be the same error. The error I’m getting is

```julia
ERROR: Zero denominator
Stacktrace:
 [1] ControlSystems.SisoRational(::ControlSystems.Poly{Float64}, ::ControlSystems.Poly{Float64}) at /Users/twan/code/discourse/v0.6/ControlSystems/src/types/sisotf.jl:7
 [2] _generic_matmatmul!(::Array{ControlSystems.SisoRational,2}, ::Char, ::Char, ::Array{ControlSystems.SisoRational,2}, ::Array{ControlSystems.SisoRational,2}) at ./linalg/matmul.jl:564
 [3] generic_matmatmul!(::Array{ControlSystems.SisoRational,2}, ::Char, ::Char, ::Array{ControlSystems.SisoRational,2}, ::Array{ControlSystems.SisoRational,2}) at ./linalg/matmul.jl:483
 [4] *(::ControlSystems.TransferFunction{ControlSystems.SisoRational}, ::ControlSystems.TransferFunction{ControlSystems.SisoRational}) at /Users/twan/code/discourse/v0.6/ControlSystems/src/types/transferfunction.jl:440
 [5] /(::ControlSystems.TransferFunction{ControlSystems.SisoRational}, ::ControlSystems.TransferFunction{ControlSystems.SisoRational}) at /Users/twan/code/discourse/v0.6/ControlSystems/src/types/transferfunction.jl:458

```

(it can be a good idea to include the full stack trace in your questions as well to improve your chances of getting quick feedback)

The issue seems to be related to  
[https://github.com/JuliaLang/julia/blob/d386e40c17d43b79fc89d3e579fc04547241787c/base/linalg/matmul.jl#L564](https://github.com/JuliaLang/julia/blob/d386e40c17d43b79fc89d3e579fc04547241787c/base/linalg/matmul.jl#L564)  
in combination with [https://github.com/JuliaControl/ControlSystems.jl/blob/5f41dcebe12b5bf20f0c88240d5d23aeb8ea6ca6/src/types/sisotf.jl#L7](https://github.com/JuliaControl/ControlSystems.jl/blob/5f41dcebe12b5bf20f0c88240d5d23aeb8ea6ca6/src/types/sisotf.jl#L7).

As a workaround, since you’re dealing with a SISO system, you could avoid the matrix multiplication code altogether and do something like

```julia
t1 = num; t2 = 1 / den
TransferFunction(fill(t1.matrix[1] * t2.matrix[1], 1, 1), t1.Ts, t2.inputnames, t1.outputnames)

```

(following [https://github.com/JuliaControl/ControlSystems.jl/blob/5f41dcebe12b5bf20f0c88240d5d23aeb8ea6ca6/src/types/transferfunction.jl#L441](https://github.com/JuliaControl/ControlSystems.jl/blob/5f41dcebe12b5bf20f0c88240d5d23aeb8ea6ca6/src/types/transferfunction.jl#L441))

By the way, I’m not super comfortable with the electrical domain so I’m not sure what `num` and `den` represent, but from a pure controls standpoint, the denominator of `num` is exactly the same as the denominator of `den`, so you’ve got pole-zero cancellation, which may be an issue. Also the denominator of the actual result is actually pretty much zero (although not exactly zero).

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 3, 2018, 5:56am UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/6 "2018-02-03T05:56:55Z")

</div>

I think it boils down to numerical issues with polynomial calculations due to the very small numerical values. Two potential solutions

- Convert your systems to statespace form, for which calculations go through.
- Checkout branch `autodiff` (PR soon to be merged at [https://github.com/JuliaControl/ControlSystems.jl/pull/118](https://github.com/JuliaControl/ControlSystems.jl/pull/118) ) which allows you to have transfer functions of `BigFloat` type, for which calculations also go through. `minreal` wont work on systems with BigFloats though which is an issue.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 4, 2018, 7:49am UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/7 "2018-02-04T07:49:51Z")

</div>

You can of course construct the transfer function from the polynomials yourself since you know there will be cancellation

```julia
tot_load_imp = tf(num.matrix[1].num.a, den.matrix[1].num.a)

```

Or you can simplify and balance your transfer functions before performing the computation

```julia
julia> den = (bulk_imp + cer_load_imp) |> minreal
TransferFunction:
0.000395831163420477s^2 + 922.6700306779511s + 659547.980638204
---------------------------------------------------------------
                   s^2 + 3719.8506107994713s

Continuous-time transfer function model

julia> num = (bulk_imp * cer_load_imp) |> minreal
TransferFunction:
0.28268226450153416s + 131909.59612764078
-----------------------------------------
        s^2 + 3719.8506107994713s

Continuous-time transfer function model

julia> tot_load_imp = num/den |> minreal
TransferFunction:
    714.1485830948866s + 3.332471223027936e8
------------------------------------------------
s^2 + 2.330968645078185e6s + 1.666235611513967e9

Continuous-time transfer function model

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 4, 2018, 8:07am UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/8 "2018-02-04T08:07:59Z")

</div>

Issue created  
[https://github.com/JuliaControl/ControlSystems.jl/issues/123](https://github.com/JuliaControl/ControlSystems.jl/issues/123)

---

<div class="post-metadata">

**Author:** ![bo-lo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bo-lo/32/3317_2.png) [@bo-lo](https://discourse.julialang.org/u/bo-lo)\
**Post date:** [February 4, 2018, 7:55pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/9 "2018-02-04T19:55:11Z")

</div>

Hi Twan,

Thanks for your email!  
In my code t1 is a s domain function.  
What does t1.matrix[1] mean?

thanks,  
Venkat.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [February 4, 2018, 8:02pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/10 "2018-02-04T20:02:56Z")

</div>

This is what `t1` (a `TransferFunction`) looks like underneath:

```julia
julia> dump(t1)
ControlSystems.TransferFunction{ControlSystems.SisoRational}
  matrix: Array{ControlSystems.SisoRational}((1, 1))
    1: ControlSystems.SisoRational
      num: ControlSystems.Poly{Float64}
        a: Array{Float64}((3,)) [5.76038e-10, 0.0002688, 0.0]
        nzfirst: Int64 1
      den: ControlSystems.Poly{Float64}
        a: Array{Float64}((4,)) [2.03776e-9, 7.58016e-6, 0.0, 0.0]
        nzfirst: Int64 1
  Ts: Float64 0.0
  nu: Int64 1
  ny: Int64 1
  inputnames: Array{String}((1,))
    1: String ""
  outputnames: Array{String}((1,))
    1: String ""

```

so `TransferFunction` is really more of a transfer function matrix, and `t1`’s matrix of `SisoRational`s is 1x1.

---

<div class="post-metadata">

**Author:** ![bo-lo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bo-lo/32/3317_2.png) [@bo-lo](https://discourse.julialang.org/u/bo-lo)\
**Post date:** [February 4, 2018, 8:30pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/11 "2018-02-04T20:30:49Z")

</div>

First of all, thanks to everybody who is responding. Amazing.

```julia
tot_load_imp = tf(num.matrix[1].num.a, den.matrix[1].num.a)

```

what does .matrix[1].num.a do? where can I find documentation on this?

---

<div class="post-metadata">

**Author:** ![bo-lo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bo-lo/32/3317_2.png) [@bo-lo](https://discourse.julialang.org/u/bo-lo)\
**Post date:** [February 4, 2018, 8:35pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/12 "2018-02-04T20:35:32Z")

</div>

Finally, this code is giving me reasonable result (0.5.2):

```julia
using ControlSystems

w = logspace(1.0, 7.5, 1000)

rl = big"0.2"

L = big"150.0e-9"
dcr = big"2.0e-3" # actual dcr is 0.2 mohm. I am adding the scaled top/bottom fet resistances.

nbulk = 12
cbulk = big"470.0e-6"
esr_bulk = big"4.5e-3"

ncer = 96
ccer = big"14.0e-6" #22uF derated for 1V DC.
esr_cer = big"2.0e-3"

nph = 1

s = tf("s")

l_imp = dcr + L*s
bulk_imp = (esr_bulk/nbulk) + 1.0/(cbulk*nbulk*s)
cer_imp = (esr_cer/ncer) + 1.0/(ccer*ncer*s)
cer_load_imp = cer_imp*rl/(cer_imp + rl) #ceramic + load resistance.
den = (bulk_imp + cer_load_imp) |> minreal
num = (bulk_imp * cer_load_imp) |> minreal

#tot_load_imp = tf(num.matrix[1].num.a, den.matrix[1].num.a)
tot_load_imp = num/den |> minreal

tot_load_tf = tot_load_imp/(l_imp + tot_load_imp) |> minreal

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 4, 2018, 8:37pm UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/13 "2018-02-04T20:37:50Z")

</div>

That line makes use of the internal representation of the polynomials in the numerator and denominator of individual SISO transfer functions. The field a is the vector of coefficients in the polynomial. Matrix[1] extracts the first and only SISO transfer function (all transfer functions are assumed to be MIMO in the toolboc). The internals of the package are currently undocumented.

---

<div class="post-metadata">

**Author:** ![bo-lo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bo-lo/32/3317_2.png) [@bo-lo](https://discourse.julialang.org/u/bo-lo)\
**Post date:** [February 6, 2018, 1:21am UTC](https://discourse.julialang.org/t/trouble-with-controlsystems-jl/8783/14 "2018-02-06T01:21:01Z")

</div>

I continued with my experiment and I have now run into issues with bodeplot. Below is the description:

```julia
using ControlSystems

w = logspace(1.0, 7.5, 1000)

Vout = big"1.0"
Iout = big"5.0"

rl = Vout/Iout

L = big"150.0e-9"
dcr = big"2.0e-3" # actual dcr is 0.2 mohm. I am adding the scaled top/bottom fet resistances.

nbulk = 12
cbulk = big"470.0e-6"
esr_bulk = big"4.5e-3"

ncer = 96
ccer = big"14.0e-6" #22uF derated for 1V DC.
esr_cer = big"2.0e-3"

nph = 1

s = tf("s")

l_imp = dcr + L*s
bulk_imp = (esr_bulk/nbulk) + 1.0/(cbulk*nbulk*s)
cer_imp = (esr_cer/ncer) + 1.0/(ccer*ncer*s)
cer_load_imp = cer_imp*rl/(cer_imp + rl) #ceramic + load resistance.
den = (bulk_imp + cer_load_imp) |> minreal
num = (bulk_imp * cer_load_imp) |> minreal

#tot_load_imp = tf(num.matrix[1].num.a, den.matrix[1].num.a)
tot_load_imp = num/den |> minreal

tot_load_tf = tot_load_imp/(l_imp + tot_load_imp) |> minreal

```

This gives me a s funtion.

```julia
bodeplot(tot_load_tf, w, title="total load tf")

```

This bode plot looked good.

```julia
zpkdata(tot_load_tf)

```

I got the poles, zeroes, gain as vectors. So far so good.  
I then visually inspected the poles and zeros and removed a matching pole/zero combination.

```julia
redu_load_tf = zpk([-4.66636e5+0.0im], [-7846.57-30034.0im,-7846.57+30034.0im,-2.32861e6-0.0im], 4.76099e9) |> minreal

```

Got the reduced s function.

```julia
bodeplot(redu_load_tf, w, title="reduced tf")

```

No bodeplot!  
I get:  
MethodError: no method matching \*(::ControlSystems.##5#7{Complex{Float64}}, ::ControlSystems.Poly{Float64})  
Closest candidates are:  
\*(::Any, ::Any, ::Any, ::Any…) at operators.jl:138  
\*(::PyCall.PyObject, ::Any) at /home/venkat/.julia/v0.5/PyCall/src/pyoperators.jl:11  
\*(::Number, ::ControlSystems.Poly{T\<:Number}) at /home/venkat/.julia/v0.5/ControlSystems/src/types/polys.jl:82  
…

in mapfoldl\_impl(::Base.#identity, ::Base.#\*, ::Function, ::Array{ControlSystems.Poly{Float64},1}, ::Int64) at ./reduce.jl:46  
in evalfr(::ControlSystems.SisoZpk, ::Complex{Float64}) at /home/venkat/.julia/v0.5/ControlSystems/src/types/sisozpk.jl:99  
in \_collect(::Array{ControlSystems.SisoZpk,2}, ::Base.Generator{Array{ControlSystems.SisoZpk,2},ControlSystems.##183#184{Complex{Float64}}}, ::Base.EltypeUnknown, ::Base.HasShape) at ./array.jl:320  
in freqresp(::ControlSystems.TransferFunction{ControlSystems.SisoZpk}, ::Array{Float64,1}) at /home/venkat/.julia/v0.5/ControlSystems/src/freqresp.jl:24  
in bode at /home/venkat/.julia/v0.5/ControlSystems/src/freqresp.jl:113 [inlined]  
in #bodeplot#208(::Bool, ::Array{Any,1}, ::Function, ::Array{ControlSystems.TransferFunction{ControlSystems.SisoZpk},1}, ::Array{Float64,1}) at /home/venkat/.julia/v0.5/ControlSystems/src/plotting.jl:141  
in (::ControlSystems.#kw##bodeplot)(::Array{Any,1}, ::ControlSystems.#bodeplot, ::Array{ControlSystems.TransferFunction{ControlSystems.SisoZpk},1}, ::Array{Float64,1}) at ./:0  
in #bodeplot#211(::Bool, ::Array{Any,1}, ::Function, ::ControlSystems.TransferFunction{ControlSystems.SisoZpk}, ::Array{Float64,1}, ::Vararg{Array{Float64,1},N}) at /home/venkat/.julia/v0.5/ControlSystems/src/plotting.jl:164  
in (::ControlSystems.#kw##bodeplot)(::Array{Any,1}, ::ControlSystems.#bodeplot, ::ControlSystems.TransferFunction{ControlSystems.SisoZpk}, ::Array{Float64,1}) at ./:0  
in include\_string(::String, ::String) at ./loading.jl:441