Class: RPowerFlow::Solver::FullNewton

Inherits:
Object
  • Object
show all
Defined in:
lib/rpowerflow/solvers/full_newton.rb

Class Method Summary collapse

Class Method Details

.solve(system, max_iterations, tolerance) ⇒ Object

Currently returns the number of iterations required to solve the system provided. Updates the system in place with the calculated values.

TODO: exceptions! Throw an exception if the solver fails to converge. Both the Linalg and NArray packages throw an exception when an LU solve is attempted against a singular matrix.

TODO: logging!



43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
 
# File 'lib/rpowerflow/solvers/full_newton.rb', line 43

def self.solve(system, max_iterations, tolerance)
  size       = system.buses.length
  iterations = 0

  max_iterations.times do
    # First half of power and delta arrays are
    # Mw values, second half are Mvar values.
    #
    # TODO: reset rather than create (possible?)
    if defined?(Linalg)
      power = Linalg::DMatrix.columns [Array.new(size * 2, 0.0)]
      delta = Linalg::DMatrix.columns [Array.new(size * 2, 0.0)]
    else
      power = NVector.float(size * 2)
      delta = NVector.float(size * 2)
    end

    jacobian = system.jacobian do |i, mw, mvar|
      bus  = system.buses[i]

      power[i]        += mw
      power[i + size] += mvar

      # As the calculated power of each bus is updated
      # above, the change in power for each bus is also
      # calculated. We could do this at the end, but it
      # would require another loop. Not sure which is
      # more expensive... doing these calculations over
      # and over again unnecessarily or doing another loop.
      # Looks like a good spot for a benchmark!
      #
      # TODO: benchmark!
      delta[i]        = bus.mw - power[i]
      delta[i + size] = bus.mvar + (bus.voltage ** 2 * bus.susceptance) - power[i + size]
    end

    tolerable = true
    (0...size).each do |i|
      bus = system.buses[i]

      # In terms of power values, Mw is a given for voltage
      # regulated, non-slack buses and Mw and Mvar are givens
      # for load buses.
      #
      # If this is a voltage regulated bus and it's not
      # the slack bus, only check to see if Mw is within
      # error range. Else, it's a load bus, so check to
      # see if both Mw and Mvar are within error range.
      if bus.avr
        unless bus.slack
          tolerable = false if delta[i].abs > tolerance
        end
      else
        tolerable = false if delta[i].abs > tolerance || delta[i + size] > tolerance
      end
    end

    # If given power values are within error range,
    # update the system buses with the calculated
    # power values and end the solve. Otherwise, solve
    # the next iteration of voltage and angle values.
    if tolerable
      (0...size).each do |i|
        bus = system.buses[i]

        bus.mw   = power[i]
        bus.mvar = power[i + size]
      end

      # Break out of the iterations loop
      # since we've solved the system.
      break
    else
      if defined?(Linalg)
        x = Linalg::DMatrix.solve(jacobian, delta)
      else
        x = jacobian.lu.solve(delta)
      end

      # Update system bus voltage and angle
      # values such that they can be used for
      # the next round of power and Jacobian
      # formulations.
      (0...size).each do |i|
        bus = system.buses[i]

        bus.voltage += x[i + size]
        bus.angle   += x[i]
      end

      # Update the iteration count since we
      # had to do another linear solve.
      iterations += 1
    end
  end

  return iterations
end