MNA solver

mna.fnl

(macros
  { :circuit
    (fn [...]
      (local k (length [...]))
      (if (= (% k 5) 0)
        (do
          (local j (/ k 5)) (local res { })
          (for [i 1 j]
            (local n (* i 5))
            (tset res i
              { :type     (. [...] (- n 4))
                :name     (. [...] (- n 3))
                :pos-node (. [...] (- n 2))
                :neg-node (. [...] (- n 1))
                :value    (. [...] n) }))
          res)
        { }))

    :replace-with
      (fn [type elem value]
        { :type type
          :value value
          :name `(. elem :name)
          :pos-node `(. elem :pos-node)
          :neg-node `(. elem :neg-node) } ) })

(local node-table {})
(tset node-table :__index node-table)
(tset node-table :create
  (fn [self]
    (local inst
      { :nodes {}
        :node-count 0
        :components 0 })
    (setmetatable inst node-table)
    inst))

(tset node-table :add-to-nodes
  (fn [self node-str]
    (when (= (. self.nodes node-str) nil)
      (tset self.nodes node-str self.node-count)
      (set self.node-count (+ self.node-count 1)))
    (. self.nodes node-str)))

(fn map-nodes [circ]
  (local tbl (node-table:create))
  (tbl:add-to-nodes "0")

  (each [_ elem (ipairs circ)]
    (tset tbl :components (+ tbl.components 1))
    (tset tbl elem.type (+ (or (. tbl elem.type) 0) 1))
    (tset elem :high (tbl:add-to-nodes elem.pos-node))
    (tset elem :low  (tbl:add-to-nodes elem.neg-node)))

  tbl)

(fn matrix-create [n m]
  (local inst { :size { :width n :height m } })

  (for [i 1 n]
    (tset inst i [])
    (for [j 1 m]
      (tset (. inst i) j 0)))
  inst)

(fn matrix-get [self i j]
  (. (or (. self i) []) j))

(fn matrix-set [self i j val]
  (tset (. self i) j val))

(fn matrix-copy [A]
  (local inst { :size { :height A.size.height
                        :width  A.size.width } })

  (for [i 1 A.size.width]
    (tset inst i [])
    (for [j 1 A.size.height]
      (tset (. inst i) j (. (. A i) j))))
  inst)

(fn matrix-print [self]
  (for [i 1 self.size.width]
    (for [j 1 self.size.height]
      (io.write (matrix-get self i j) " "))
    (io.write "\n")))

(fn calculate-resistor [A b elem g2-index]
  (when (~= elem.high 0)
    (matrix-set A elem.high elem.high
      (+ (matrix-get A elem.high elem.high)
         (/ 1 elem.value))))
  (when (~= elem.low 0)
    (matrix-set A elem.low elem.low
      (+ (matrix-get A elem.low elem.low)
         (/ 1 elem.value))))
  (when (and (~= elem.high 0) (~= elem.low 0))
    (matrix-set A elem.high elem.low
      (- (matrix-get A elem.high elem.low)
         (/ 1 elem.value)))
    (matrix-set A elem.low elem.high
      (- (matrix-get A elem.low elem.high)
         (/ 1 elem.value)))))

(fn calculate-voltage [A b elem g2-index]
  (when (~= elem.high 0)
    (matrix-set A elem.high g2-index
      (+ (matrix-get A elem.high g2-index) 1))
    (matrix-set A g2-index elem.high
      (+ (matrix-get A g2-index elem.high) 1)))

  (when (~= elem.low 0)
    (matrix-set A elem.low g2-index
      (- (matrix-get A elem.low g2-index) 1))
    (matrix-set A g2-index elem.low
      (- (matrix-get A g2-index elem.low) 1)))

  (matrix-set b g2-index 1 elem.value)
  (+ g2-index 1))

(fn calculate-current [A b elem g2-index]
  (when (~= elem.high 0)
    (matrix-set b elem.high 1 (- (matrix-get b elem.high 1) elem.value)))
  (when (~= elem.low 0)
    (matrix-set b elem.low 1 (+ (matrix-get b elem.low 1) elem.value))))

(fn calculate-inductance [A b elem g2-index]
  (when (~= elem.high 0)
    (matrix-set A elem.high g2-index
      (+ (matrix-get A elem.high g2-index) 1))
    (matrix-set A g2-index elem.high
      (+ (matrix-get A g2-index elem.high) 1)))

  (when (~= elem.low 0)
    (matrix-set A elem.low g2-index
      (- (matrix-get A elem.low g2-index) 1))
    (matrix-set A g2-index elem.low
      (- (matrix-get A g2-index elem.low) 1)))

  (matrix-set b g2-index 1 0)
  (+ g2-index 1))

(fn calculate-capacitor [A b elem g2-index])

(local circuit-elems
  { :resistor  calculate-resistor
    :voltage   calculate-voltage
    :current   calculate-current
    :inductor  calculate-inductance
    :capacitor calculate-capacitor })

(require :solve)

(fn solve-dc-aux [tbl circ]
  (let [g2-count (+ tbl.voltage (or tbl.inductor 0))
        matrix-size (+ tbl.node-count g2-count -1)]

    (var A (matrix-create matrix-size matrix-size))
    (var b (matrix-create matrix-size 1))

    (var g2-index (- matrix-size g2-count -1))
    ;; generate A and b
    (each [id elem (ipairs circ)]
      (let [func (. circuit-elems elem.type)]
        (local maybe-g2-index (func A b elem g2-index))
        (when (~= maybe-g2-index nil)
          (set g2-index maybe-g2-index)
          (tset elem :current-index (- maybe-g2-index 1)))))

    (local solution (solve A b))

    (var res { :voltages {} :currents {} } )
    (each [name pin (pairs tbl.nodes)]
      (let [v (or (matrix-get solution pin 1) 0)]
        (tset res.voltages pin v)))

    (each [id elem (ipairs circ)]
      (when (~= elem.current-index nil)
        (tset res.currents id
          (matrix-get solution elem.current-index 1))))

    res))

(fn solve-dc [circ] (solve-dc-aux (map-nodes circ) circ))

(fn concat [tbl delim]
  (var res "")

  (each [id val (pairs tbl)]
    (when (~= res "") (set res (.. res ";")))
    (set res (.. res val)))

  res)

(local ac-elems [ :inductance :capacitance :ac-sources ])
(fn solve-ac [circ Δt T]
  (var vars {}) (var circ-reduced {})
  (each [_ name (ipairs ac-elems)]
    (tset vars name {}))

  (each [id elem (ipairs circ)]
    (match elem.type
      :inductor (do
        (tset circ-reduced id (replace-with :voltage elem 0))
        (tset vars.inductance id { :L elem.value :I 0 }))
      :capacitor (do
        (tset circ-reduced id (replace-with :current elem 0))
        (tset vars.capacitance id { :C elem.value :U 0 }))
      :ac-source (do
        (tset circ-reduced id
          (replace-with :voltage elem (elem.value 0)))
        (tset vars.ac-sources id elem.value))
      _ (tset circ-reduced id elem)))

  (let [tbl (map-nodes circ-reduced)]
    (var res (solve-dc circ-reduced))
    ;; init inductors current
    (each [id elem (pairs vars.inductance)]
      (set elem.I (. res.currents id)))
    ;; init capacitor voltages
    (each [id elem (pairs vars.capacitance)]
      (let [low (. circ-reduced id :low)
            high (. circ-reduced id :high)
            U-low (. res.voltages low)
            U-high (. res.voltages high)
            U (- U-high U-low)]
        (set (elem.U elem.low elem.high)
             (values U low high))))

    (for [t 0 T Δt]
      ;; recalculate as for DC
      (set res (solve-dc circ-reduced))

      ;; update emf of inductors
      (each [id elem (pairs vars.inductance)]
        (let [I (. res.currents id)
              ΔI (- I elem.I)
              ε (* (- elem.L) (/ ΔI Δt))]
          (set elem.I I) (tset (. circ-reduced id) :value ε)))

      ;; update current of capacitors
      (each [id elem (pairs vars.capacitance)]
        (let [U-low (. res.voltages elem.low)
              U-high (. res.voltages elem.high)
              U (- U-high U-low)
              ΔU (- U elem.U)
              I (* elem.C (/ ΔU Δt))]
          (set elem.U U) (tset (. circ-reduced id) :value I)))

      ;; update emf of AC sources
      (each [id ε (pairs vars.ac-sources)]
        (tset (. circ-reduced id) :value (ε t)))

      (print (.. t ";" (concat res.currents ";") ";" (concat res.voltages ";"))))

    res))

(local circ-0
  (circuit
    :voltage  :A :1 :0 10
    :resistor :B :1 :2 3
    :resistor :C :2 :3 2
    :voltage  :X :3 :0 0))

(local circ-1
  (circuit
    :voltage   :A :5  :0 2
    :voltage   :B :3  :2 0.2
    :voltage   :C :7  :6 2
    :resistor  :D :1  :5 1.5
    :resistor  :E :1  :12 1
    :resistor  :F :5  :2 50
    :resistor  :G :5  :6 0.1
    :resistor  :H :2  :6 1.5
    :resistor  :I :3  :4 0.2
    :resistor  :J :7  :0 1e3
    :resistor  :K :4  :0 10
    :current   :L :4  :7 1e-3
    :current   :M :0  :6 1e-3
    :capacitor :N :7  :0 0.1
    :capacitor :X :2  :0 0.2
    :inductor  :Y :12 :2 0.1))

(fn sinusoidal-ac-source [εₘ ω]
  (fn [t] (* εₘ (math.sin (* ω t)))))

(local circ-2
  (circuit
    :ac-source :A :1 :0 (sinusoidal-ac-source 50 (* 2 math.pi 1))
    :resistor  :B :1 :2 100
    :capacitor :D :2 :0 0.0001
    :inductor  :C :2 :0 1))

(local solution (solve-ac circ-2 0.1 60))

solve.lua

-- from https://rosettacode.org/wiki/Gaussian_elimination#Python
-- crushes “a” and “b”
function solve(a, b)
    local n = a.size.width
    local p = b.size.height

    for i = 1, n do
        local k = i
        for j = i + 1, n do
            if math.abs(a[j][i]) > math.abs(a[k][i]) then
                k = j
            end
        end

        if k ~= i then
            a[i], a[k] = a[k], a[i]
            b[i], b[k] = b[k], b[i]
        end

        for j = i + 1, n do
            local t = a[j][i] / a[i][i]
            for k = i + 1, n do
                a[j][k] = a[j][k] - t * a[i][k]
            end

            for k = 1, p do
                b[j][k] = b[j][k] - t * b[i][k]
            end
        end
    end

    for i = n, 1, -1 do
        for j = i + 1, n do
            local t = a[i][j]
            for k = 1, p do
                b[i][k] = b[i][k] - t * b[j][k]
            end
        end

        local t = 1/a[i][i]
        for j = 1, p do
            b[i][j] = b[i][j] * t
        end
    end

    return b
end