Trouble with ControlSystems.jl

bo-lo · February 2, 2018, 5:37pm

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):

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.

carstenbauer · February 2, 2018, 11:39pm

Please have a look at: PSA: how to quote code with backticks

baggepinnen · February 3, 2018, 12:21am

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

bo-lo · February 3, 2018, 12:31am

Thanks,
I have updated my post.

tkoolen · February 3, 2018, 5:07am

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

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
in combination with 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

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)

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).

baggepinnen · February 3, 2018, 5:56am

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 ) 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.

baggepinnen · February 4, 2018, 7:49am

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

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> 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

baggepinnen · February 4, 2018, 8:07am

Issue created
https://github.com/JuliaControl/ControlSystems.jl/issues/123

bo-lo · February 4, 2018, 7:55pm

Hi Twan,

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

thanks,
Venkat.

tkoolen · February 4, 2018, 8:02pm

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

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 SisoRationals is 1x1.

bo-lo · February 4, 2018, 8:30pm

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

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?

bo-lo · February 4, 2018, 8:35pm

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

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

baggepinnen · February 4, 2018, 8:37pm

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.

bo-lo · February 6, 2018, 1:21am

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

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.

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

This bode plot looked good.

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.

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.

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