Opened 5 years ago

Last modified 2 years ago

#5302 new defect

PowerSystems models with redundant intial conditions cannot be handled by backend

Reported by: casella Owned by: lochel
Priority: blocker Milestone: 2.0.0
Component: Backend Version:
Keywords: Cc: Karim.Abdelhak, AnHeuermann, rfranke

Description

The following models in the PowerSystems library

PowerSystems.Examples.AC1ph_DC.Inverters.InverterToGrid
PowerSystems.Examples.AC3ph.Generation.PowerAngle
PowerSystems.Examples.AC3ph.Generation.Vsource
PowerSystems.Examples.AC3ph.Inverters.InverterToGrid
PowerSystems.Examples.AC3ph.Transformation.PhaseShifts
PowerSystems.Examples.AC3ph.Transmission.DoubleRXline
PowerSystems.Examples.AC3ph.Transmission.DoublelLine

have redundant consistent initial conditions, as diagnosed by Dymola. OMC has some heuristics to detect them and remove them, but apparently we have a problem handling this case in OMC.

Regarding the first option, I analyzed the d=optdaedump output of PowerSystems.Examples.AC3ph.Generation.Vsource. After running Pantelides' algorithm, there are 26 variables/equations:

1: grd2.term.i:VARIABLE(flow=true unit = "A" )  "Current flowing into the pin" annotation(unassignedMessage = "An electrical current cannot be uniquely calculated. 
The reason could be that 
- a ground object is missing (Modelica.Electrical.Analog.Basic.Ground) 
  to define the zero potential of the electrical circuit, or 
- a connector of an electrical component is not connected.") type: Real  
2: grd1.term.i:VARIABLE(flow=true unit = "A" )  "Current flowing into the pin" annotation(unassignedMessage = "An electrical current cannot be uniquely calculated. 
The reason could be that 
- a ground object is missing (Modelica.Electrical.Analog.Basic.Ground) 
  to define the zero potential of the electrical circuit, or 
- a connector of an electrical component is not connected.") type: Real  
3: infBus.theta:DUMMY_STATE(flow=false start = 0.0 unit = "rad" stateSelect=StateSelect.prefer protected = true )  "absolute angle" type: Real  
4: line.omega[2]:VARIABLE(flow=false unit = "rad/s" protected = true )  type: Real  [2] 
5: line.i[3]:STATE(1)(flow=false start = line.i_start[3] unit = "A" nominal = 1.0 )  type: Real  [3] 
6: line.i[2]:STATE(1)(flow=false start = line.i_start[2] unit = "A" nominal = 1.0 )  type: Real  [3] 
7: line.i[1]:STATE(1)(flow=false start = line.i_start[1] unit = "A" nominal = 1.0 )  type: Real  [3] 
8: line.v[2]:VARIABLE(flow=false start = line.v_start[2] unit = "V" nominal = 1000.0 )  type: Real  [3] 
9: line.v[1]:VARIABLE(flow=false start = line.v_start[1] unit = "V" nominal = 1000.0 )  type: Real  [3] 
10: line.term_n.v[2]:VARIABLE(flow=false unit = "V" nominal = 1000.0 )  "voltage vector" type: Real  [3] 
11: line.term_n.v[1]:VARIABLE(flow=false unit = "V" nominal = 1000.0 )  "voltage vector" type: Real  [3] 
12: sensor.p[3]:VARIABLE(flow=false )  "{active, reactive, DC} power, term_p to term_n" annotation(Placement(transformation(origin = {0, 100}, extent = {{-10, -10}, {10, 10}}, rotation = 90))) type: Real  [3] 
13: sensor.p[2]:VARIABLE(flow=false )  "{active, reactive, DC} power, term_p to term_n" annotation(Placement(transformation(origin = {0, 100}, extent = {{-10, -10}, {10, 10}}, rotation = 90))) type: Real  [3] 
14: sensor.p[1]:VARIABLE(flow=false )  "{active, reactive, DC} power, term_p to term_n" annotation(Placement(transformation(origin = {0, 100}, extent = {{-10, -10}, {10, 10}}, rotation = 90))) type: Real  [3] 
15: Vsource.omega[2]:VARIABLE(flow=false unit = "rad/s" protected = true )  type: Real  [2] 
16: Vsource.omega[1]:VARIABLE(flow=false unit = "rad/s" protected = true )  type: Real  [2] 
17: Vsource.i[3]:DUMMY_STATE(flow=false start = 0.0 unit = "A" nominal = 1.0 protected = true )  type: Real  [3] 
18: Vsource.i[2]:DUMMY_STATE(flow=false start = 0.0 unit = "A" nominal = 1.0 protected = true )  type: Real  [3] 
19: Vsource.i[1]:DUMMY_STATE(flow=false start = 0.0 unit = "A" nominal = 1.0 protected = true )  type: Real  [3] 
20: Vsource.v[3]:VARIABLE(flow=false start = line.v_start[3] unit = "V" nominal = 1000.0 protected = true )  type: Real  [3] 
21: Vsource.v[2]:VARIABLE(flow=false start = 0.0 unit = "V" nominal = 1000.0 protected = true )  type: Real  [3] 
22: Vsource.v[1]:VARIABLE(flow=false start = 20000.0 unit = "V" nominal = 1000.0 protected = true )  type: Real  [3] 
23: $DER.Vsource.i[1]:DUMMY_DER(flow=false fixed = false )  type: Real  [3] 
24: $DER.Vsource.i[2]:DUMMY_DER(flow=false fixed = false )  type: Real  [3] 
25: $DER.Vsource.i[3]:DUMMY_DER(flow=false fixed = false )  type: Real  [3] 
26: $DER.infBus.theta:DUMMY_DER(flow=false fixed = false )  "absolute angle" type: Real  
 
 
Equations (26, 26) 
======================================== 
1/1 (1): $DER.infBus.theta = 314.1592653589793   [dynamic |0|0|0|0|] 
2/2 (1): line.omega[2] = Vsource.omega[2]   [dynamic |0|0|0|0|] 
3/3 (1): Vsource.L * $DER.Vsource.i[1] + Vsource.R * Vsource.i[1] - Vsource.i[2] * Vsource.omega[2] * Vsource.L = Vsource.V * cos(Vsource.alpha0) - Vsource.v[1]   [dynamic |0|0|0|0|] 
4/4 (1): Vsource.L * $DER.Vsource.i[2] + Vsource.i[1] * Vsource.omega[2] * Vsource.L + Vsource.R * Vsource.i[2] = Vsource.V * sin(Vsource.alpha0) - Vsource.v[2]   [dynamic |0|0|0|0|] 
5/5 (1): Vsource.L0 * $DER.Vsource.i[3] + Vsource.R * Vsource.i[3] = -Vsource.v[3]   [dynamic |0|0|0|0|] 
6/6 (1): (-1.732050807568877) * line.i[3] - grd1.term.i = 0.0   [dynamic |0|0|0|0|] 
7/7 (1): sensor.p[1] = Vsource.v[1] * line.i[1] + Vsource.v[2] * line.i[2]   [dynamic |0|0|0|0|] 
8/8 (1): sensor.p[2] = Vsource.v[2] * line.i[1] - Vsource.v[1] * line.i[2]   [dynamic |0|0|0|0|] 
9/9 (1): sensor.p[3] = Vsource.v[3] * line.i[3]   [dynamic |0|0|0|0|] 
10/10 (1): line.omega[2] = $DER.infBus.theta   [dynamic |0|0|0|0|] 
11/11 (1): line.v[1] = Vsource.v[1] - line.term_n.v[1]   [dynamic |0|0|0|0|] 
12/12 (1): line.v[2] = Vsource.v[2] - line.term_n.v[2]   [dynamic |0|0|0|0|] 
13/13 (1): line.L * der(line.i[1]) - line.i[2] * line.omega[2] * line.L + line.R * line.i[1] = line.v[1]   [dynamic |0|0|0|0|] 
14/14 (1): line.L * der(line.i[2]) + line.i[1] * line.omega[2] * line.L + line.R * line.i[2] = line.v[2]   [dynamic |0|0|0|0|] 
15/15 (1): line.L0 * der(line.i[3]) + line.R * line.i[3] = Vsource.v[3]   [dynamic |0|0|0|0|] 
16/16 (1): line.term_n.v[1] = infBus.V * cos(infBus.alpha0)   [dynamic |0|0|0|0|] 
17/17 (1): line.term_n.v[2] = infBus.V * sin(infBus.alpha0)   [dynamic |0|0|0|0|] 
18/18 (1): 1.732050807568877 * line.i[3] - grd2.term.i = 0.0   [dynamic |0|0|0|0|] 
19/19 (1): Vsource.omega[1] = line.omega[1]   [dynamic |0|0|0|0|] 
20/20 (1): $DER.Vsource.i[3] = der(line.i[3])   [dynamic |0|0|0|0|] 
21/21 (1): $DER.Vsource.i[2] = der(line.i[2])   [dynamic |0|0|0|0|] 
22/22 (1): $DER.Vsource.i[1] = der(line.i[1])   [dynamic |0|0|0|0|] 
23/23 (1): Vsource.i[1] = line.i[1]   [dynamic |0|0|0|0|] 
24/24 (1): Vsource.i[2] = line.i[2]   [dynamic |0|0|0|0|] 
25/25 (1): Vsource.i[3] = line.i[3]   [dynamic |0|0|0|0|] 
26/26 (1): infBus.theta = 314.1592653589793 * time   [dynamic |0|0|0|0|] 

There are 4 initial equations:

Initial Equations (4, 8) 
======================================== 
1/1 (1): sqrt(Vsource.v[1] ^ 2.0 + Vsource.v[2] ^ 2.0) = 20000.0 * Vsource.v_start   [initial |0|0|0|0|] 
2/2 (1): PowerSystems.AC3ph.Sources.VsourceRX.atan2(Vsource.v[2], Vsource.v[1]) = Vsource.alpha_start   [initial |0|0|0|0|] 
3/3 (3): der(Vsource.i) = {(-Vsource.i[2]) * Vsource.omega[1], Vsource.i[1] * Vsource.omega[1], 0.0}   [dynamic |0|0|0|0|] 
4/6 (3): der(line.i) = {(-line.i[2]) * line.omega[1], line.i[1] * line.omega[1], 0.0}   [dynamic |0|0|0|0|] 

The last two equations are three-dimensional array equations, so we have a total of eight scalar initial equations.

When the initial system of equations is finally created, it is reported as:

######################################## 
created initial system 
######################################## 
 
 
unknown partition 
======================================== 
 
Variables (44) 
======================================== 
1: Vsource.v0:VARIABLE(flow=false start = 1.0 unit = "V/V" fixed = false protected = true final = true )  "voltage behind reactance" type: Real  
2: Vsource.V:VARIABLE(flow=false start = 20000.0 unit = "V" nominal = 1000.0 protected = true final = true )  = 20000.0 * Vsource.v0  type: Real  
3: Vsource.alpha0:VARIABLE(flow=false start = 0.0 unit = "rad" fixed = false protected = true final = true )  "phase angle of voltage b r" type: Real  
4: grd2.term.v:VARIABLE(flow=false unit = "V" )  = 0.0  "Potential at the pin" annotation(unassignedMessage = "An electrical potential cannot be uniquely calculated. 
The reason could be that 
- a ground object is missing (Modelica.Electrical.Analog.Basic.Ground) 
  to define the zero potential of the electrical circuit, or 
- a connector of an electrical component is not connected.") type: Real  
5: grd1.term.v:VARIABLE(flow=false unit = "V" )  = 0.0  "Potential at the pin" annotation(unassignedMessage = "An electrical potential cannot be uniquely calculated. 
The reason could be that 
- a ground object is missing (Modelica.Electrical.Analog.Basic.Ground) 
  to define the zero potential of the electrical circuit, or 
- a connector of an electrical component is not connected.") type: Real  
6: system.omega:VARIABLE(flow=false start = 314.1592653589793 unit = "rad/s" )  = 314.1592653589793  type: Real  
7: infBus.V:VARIABLE(flow=false unit = "V" nominal = 1000.0 protected = true )  = 20000.0 * infBus.v0  type: Real  
8: line.term_n.v[3]:VARIABLE(flow=false unit = "V" nominal = 1000.0 )  = 0.0  "voltage vector" type: Real  [3] 
9: line.term_p.theta[1]:DUMMY_STATE(flow=false )  = 0.0  "optional vector of phase angles" type: Real  [2] 
10: line.omega[1]:VARIABLE(flow=false unit = "rad/s" protected = true )  = 0.0  type: Real  [2] 
11: system.receiveFreq.w_h:VARIABLE(flow=false )  = 0.0  "Dummy potential-variable to balance flow-variable w_H" type: Real  
12: system.receiveFreq.h:VARIABLE(flow=false )  = 0.0  "Dummy potential-variable to balance flow-variable H" type: Real  
13: system.receiveFreq.w_H:VARIABLE(flow=true unit = "rad" )  = 0.0  "angular velocity, inertia-weighted" type: Real  
14: system.receiveFreq.H:VARIABLE(flow=true unit = "s" )  = 0.0  "inertia constant" type: Real  
15: grd2.term.i:VARIABLE(flow=true unit = "A" )  "Current flowing into the pin" annotation(unassignedMessage = "An electrical current cannot be uniquely calculated. 
The reason could be that 
- a ground object is missing (Modelica.Electrical.Analog.Basic.Ground) 
  to define the zero potential of the electrical circuit, or 
- a connector of an electrical component is not connected.") type: Real  
16: grd1.term.i:VARIABLE(flow=true unit = "A" )  "Current flowing into the pin" annotation(unassignedMessage = "An electrical current cannot be uniquely calculated. 
The reason could be that 
- a ground object is missing (Modelica.Electrical.Analog.Basic.Ground) 
  to define the zero potential of the electrical circuit, or 
- a connector of an electrical component is not connected.") type: Real  
17: infBus.theta:DUMMY_STATE(flow=false start = 0.0 unit = "rad" stateSelect=StateSelect.prefer protected = true )  "absolute angle" type: Real  
18: line.omega[2]:VARIABLE(flow=false unit = "rad/s" protected = true )  type: Real  [2] 
19: $DER.line.i[3]:VARIABLE(flow=false start = line.i_start[3] unit = "A" nominal = 1.0 )  type: Real  [3] 
20: line.i[3]:VARIABLE(flow=false start = line.i_start[3] unit = "A" nominal = 1.0 )  type: Real  [3] 
21: $DER.line.i[2]:VARIABLE(flow=false start = line.i_start[2] unit = "A" nominal = 1.0 )  type: Real  [3] 
22: line.i[2]:VARIABLE(flow=false start = line.i_start[2] unit = "A" nominal = 1.0 )  type: Real  [3] 
23: $DER.line.i[1]:VARIABLE(flow=false start = line.i_start[1] unit = "A" nominal = 1.0 )  type: Real  [3] 
24: line.i[1]:VARIABLE(flow=false start = line.i_start[1] unit = "A" nominal = 1.0 )  type: Real  [3] 
25: line.v[2]:VARIABLE(flow=false start = line.v_start[2] unit = "V" nominal = 1000.0 )  type: Real  [3] 
26: line.v[1]:VARIABLE(flow=false start = line.v_start[1] unit = "V" nominal = 1000.0 )  type: Real  [3] 
27: line.term_n.v[2]:VARIABLE(flow=false unit = "V" nominal = 1000.0 )  "voltage vector" type: Real  [3] 
28: line.term_n.v[1]:VARIABLE(flow=false unit = "V" nominal = 1000.0 )  "voltage vector" type: Real  [3] 
29: sensor.p[3]:VARIABLE(flow=false )  "{active, reactive, DC} power, term_p to term_n" annotation(Placement(transformation(origin = {0, 100}, extent = {{-10, -10}, {10, 10}}, rotation = 90))) type: Real  [3] 
30: sensor.p[2]:VARIABLE(flow=false )  "{active, reactive, DC} power, term_p to term_n" annotation(Placement(transformation(origin = {0, 100}, extent = {{-10, -10}, {10, 10}}, rotation = 90))) type: Real  [3] 
31: sensor.p[1]:VARIABLE(flow=false )  "{active, reactive, DC} power, term_p to term_n" annotation(Placement(transformation(origin = {0, 100}, extent = {{-10, -10}, {10, 10}}, rotation = 90))) type: Real  [3] 
32: Vsource.omega[2]:VARIABLE(flow=false unit = "rad/s" protected = true )  type: Real  [2] 
33: Vsource.omega[1]:VARIABLE(flow=false unit = "rad/s" protected = true )  type: Real  [2] 
34: Vsource.i[3]:DUMMY_STATE(flow=false start = 0.0 unit = "A" nominal = 1.0 protected = true )  type: Real  [3] 
35: Vsource.i[2]:DUMMY_STATE(flow=false start = 0.0 unit = "A" nominal = 1.0 protected = true )  type: Real  [3] 
36: Vsource.i[1]:DUMMY_STATE(flow=false start = 0.0 unit = "A" nominal = 1.0 protected = true )  type: Real  [3] 
37: Vsource.v[3]:VARIABLE(flow=false start = line.v_start[3] unit = "V" nominal = 1000.0 protected = true )  type: Real  [3] 
38: Vsource.v[2]:VARIABLE(flow=false start = 0.0 unit = "V" nominal = 1000.0 protected = true )  type: Real  [3] 
39: Vsource.v[1]:VARIABLE(flow=false start = 20000.0 unit = "V" nominal = 1000.0 protected = true )  type: Real  [3] 
40: $DER.Vsource.i[1]:DUMMY_DER(flow=false fixed = false )  type: Real  [3] 
41: $DER.Vsource.i[2]:DUMMY_DER(flow=false fixed = false )  type: Real  [3] 
42: $DER.Vsource.i[3]:DUMMY_DER(flow=false fixed = false )  type: Real  [3] 
43: $DER.infBus.theta:DUMMY_DER(flow=false fixed = false )  "absolute angle" type: Real  
44: system.initime:DISCRETE(flow=false unit = "s" fixed = false )  type: Real  
 
 
Equations (43, 47) 
======================================== 
1/1 (1): sqrt(Vsource.v[1] ^ 2.0 + Vsource.v[2] ^ 2.0) = 20000.0 * Vsource.v_start   [initial |0|0|0|0|] 
2/2 (1): PowerSystems.AC3ph.Sources.VsourceRX.atan2(Vsource.v[2], Vsource.v[1]) = Vsource.alpha_start   [initial |0|0|0|0|] 
3/3 (3): {$DER.Vsource.i[1], $DER.Vsource.i[2], $DER.Vsource.i[3]} = {(-Vsource.i[2]) * Vsource.omega[1], Vsource.i[1] * Vsource.omega[1], 0.0}   [dynamic |0|0|0|0|] 
4/6 (3): {$DER.line.i[1], $DER.line.i[2], $DER.line.i[3]} = {(-line.i[2]) * line.omega[1], line.i[1] * line.omega[1], 0.0}   [dynamic |0|0|0|0|] 
5/9 (1): infBus.theta = 314.1592653589793 * time   [dynamic |0|0|0|0|] 
6/10 (1): Vsource.i[3] = line.i[3]   [dynamic |0|0|0|0|] 
7/11 (1): Vsource.i[2] = line.i[2]   [dynamic |0|0|0|0|] 
8/12 (1): Vsource.i[1] = line.i[1]   [dynamic |0|0|0|0|] 
9/13 (1): $DER.Vsource.i[1] = $DER.line.i[1]   [dynamic |0|0|0|0|] 
10/14 (1): $DER.Vsource.i[2] = $DER.line.i[2]   [dynamic |0|0|0|0|] 
11/15 (1): $DER.Vsource.i[3] = $DER.line.i[3]   [dynamic |0|0|0|0|] 
12/16 (1): Vsource.omega[1] = line.omega[1]   [dynamic |0|0|0|0|] 
13/17 (1): 1.732050807568877 * line.i[3] - grd2.term.i = 0.0   [dynamic |0|0|0|0|] 
14/18 (1): line.term_n.v[2] = infBus.V * sin(infBus.alpha0)   [dynamic |0|0|0|0|] 
15/19 (1): line.term_n.v[1] = infBus.V * cos(infBus.alpha0)   [dynamic |0|0|0|0|] 
16/20 (1): line.L0 * $DER.line.i[3] + line.R * line.i[3] = Vsource.v[3]   [dynamic |0|0|0|0|] 
17/21 (1): line.L * $DER.line.i[2] + line.i[1] * line.omega[2] * line.L + line.R * line.i[2] = line.v[2]   [dynamic |0|0|0|0|] 
18/22 (1): line.L * $DER.line.i[1] - line.i[2] * line.omega[2] * line.L + line.R * line.i[1] = line.v[1]   [dynamic |0|0|0|0|] 
19/23 (1): line.v[2] = Vsource.v[2] - line.term_n.v[2]   [dynamic |0|0|0|0|] 
20/24 (1): line.v[1] = Vsource.v[1] - line.term_n.v[1]   [dynamic |0|0|0|0|] 
21/25 (1): line.omega[2] = $DER.infBus.theta   [dynamic |0|0|0|0|] 
22/26 (1): sensor.p[3] = Vsource.v[3] * line.i[3]   [dynamic |0|0|0|0|] 
23/27 (1): sensor.p[2] = Vsource.v[2] * line.i[1] - Vsource.v[1] * line.i[2]   [dynamic |0|0|0|0|] 
24/28 (1): sensor.p[1] = Vsource.v[1] * line.i[1] + Vsource.v[2] * line.i[2]   [dynamic |0|0|0|0|] 
25/29 (1): (-1.732050807568877) * line.i[3] - grd1.term.i = 0.0   [dynamic |0|0|0|0|] 
26/30 (1): Vsource.L0 * $DER.Vsource.i[3] + Vsource.R * Vsource.i[3] = -Vsource.v[3]   [dynamic |0|0|0|0|] 
27/31 (1): Vsource.L * $DER.Vsource.i[2] + Vsource.i[1] * Vsource.omega[2] * Vsource.L + Vsource.R * Vsource.i[2] = Vsource.V * sin(Vsource.alpha0) - Vsource.v[2]   [dynamic |0|0|0|0|] 
28/32 (1): Vsource.L * $DER.Vsource.i[1] + Vsource.R * Vsource.i[1] - Vsource.i[2] * Vsource.omega[2] * Vsource.L = Vsource.V * cos(Vsource.alpha0) - Vsource.v[1]   [dynamic |0|0|0|0|] 
29/33 (1): line.omega[2] = Vsource.omega[2]   [dynamic |0|0|0|0|] 
30/34 (1): $DER.infBus.theta = 314.1592653589793   [dynamic |0|0|0|0|] 
31/35 (1): system.initime = time   [dynamic |0|0|0|0|] 
32/36 (1): Vsource.V = 20000.0 * Vsource.v0   [initial |0|0|0|0|] 
33/37 (1): grd2.term.v = 0.0   [initial |0|0|0|0|] 
34/38 (1): grd1.term.v = 0.0   [initial |0|0|0|0|] 
35/39 (1): system.omega = 314.1592653589793   [initial |0|0|0|0|] 
36/40 (1): infBus.V = 20000.0 * infBus.v0   [initial |0|0|0|0|] 
37/41 (1): line.term_n.v[3] = 0.0   [initial |0|0|0|0|] 
38/42 (1): line.term_p.theta[1] = 0.0   [initial |0|0|0|0|] 
39/43 (1): line.omega[1] = 0.0   [initial |0|0|0|0|] 
40/44 (1): system.receiveFreq.w_h = 0.0   [initial |0|0|0|0|] 
41/45 (1): system.receiveFreq.h = 0.0   [initial |0|0|0|0|] 
42/46 (1): system.receiveFreq.w_H = 0.0   [initial |0|0|0|0|] 
43/47 (1): system.receiveFreq.H = 0.0   [initial |0|0|0|0|] 

Here, I cannot follow exactly what the back-end is doing. On one hand, additional equations and variables have been added on top of the 26+8 that I expected at the end of Pantelides. On the other hand, it seems that the two array initial equations have disappeared. In fact, three of them are redundant, because of index reduction, but the other three are not. Eventually the initial system is unbalanced, but there are a number of fishy things going on here, and I'm not sure what is the root cause of the issue with these models.

I need someone proficient with the backend to continue the analysis. I suspect that the three-dimensional array initial equations are somehow lost, but I'm not sure if I am interpreting the optdaedump output correctly.

Change History (3)

comment:1 Changed 5 years ago by casella

  • Cc rfranke added

comment:2 Changed 5 years ago by casella

  • Cc Karim.Abdelhak AnHeuermann added
  • Milestone changed from 1.14.0 to 2.0.0
  • Priority changed from high to blocker

Karim, Andreas, this should be resolved for 2.0.0

comment:3 Changed 2 years ago by casella

Most of the initially mentioned cases now work fine, thanks to the improvements to the ASCC algorithm. There are a few remaining broken models:

Probably related to #6196

Note: See TracTickets for help on using tickets.