Solving time-dependent system of equations from Modelica with ModelingToolkit - hasmetadata error

I’m trying to find \bar{v}_{dq}^+ and \bar{v}_{dq}^- in the equation below


where {v}_{dq}^+ and {v}_{dq}^- are known time-varying signals. I’ve done this in Modelica, but I’d like to try with Julia. This is the Modelica code, where xdq1=v_{dq}^+, and xdq1M=\bar{v}_{dq}^+:

model DCinitSimp
  import Modelica;
  import Modelica.Constants.pi;

  parameter Real V1 = 0.9 "Input phase 1 voltage";
  parameter Real V2 = 1 "Input phase 2 voltage";
  parameter Real V3 = 1.05 "Input phase 3 voltage";
  parameter Real f = 50 "Grid frequency";
  parameter Real Phi = 1 "Grid phase shift";

  output Real xdq1M[2], xdq2M[2];
  
  Real x1, x2, x3;
  Real theta;
  Real xdq1[2], xdq2[2];
  Real Theta[2,2],minusTheta[2,2];
  Real Eye[2,2];

equation
  Eye = 0.00001*[1, 0; 0, 1];
  theta = 2*pi*f*time;
  Theta = [cos(theta), sin(theta); -sin(theta), cos(theta)]+Eye;
  minusTheta = [cos(-theta), sin(-theta); -sin(-theta), cos(-theta)]+Eye;
  
  x1 = V1 * cos(2*pi*f*time - Phi);
  x2 = V2 * cos(2*pi*f*time - (2/3)*pi - Phi);
  x3 = V3 * cos(2*pi*f*time + (2/3)*pi - Phi);
  
  xdq1 = Theta * sqrt(2/3) * [1, -0.5, -0.5; 0, sqrt(3)/2, -sqrt(3)/2] * {x1, x2, x3};
  xdq2 = minusTheta * sqrt(2/3) * [1, -0.5, -0.5; 0, sqrt(3)/2, -sqrt(3)/2] * {x1, x2, x3};
  
  xdq1 = xdq1M +  Theta * Theta * xdq2M;
  xdq2 = xdq2M +  minusTheta * minusTheta * xdq1M;

annotation (
     preferredView = "diagram",
     experiment(StartTime = 0.0, StopTime = 0.1, Tolerance = 1e-06, Interval = 0.0002));
end DCinitSimp;

And this is the translation to Julia so far

using Symbolics, NLsolve, ModelingToolkit

V1 = 0.9
V2 = 1
V3 = 1.05
f = 50
Φ = 1
@parameters t
@variables θ Theta[1:2,1:2] minusTheta[1:2,1:2] x1 x2 x3 xdq1[1:2] xdq2[1:2] xdq1M[1:2] xdq2M[1:2] 

eqs = [
       θ ~ t,
       x1 ~ V1 * cos(2*pi*f*t- Φ),
       x2 ~ V2 * cos(2*pi*f*t- (2/3)*pi - Φ),
       x3 ~ V3 * cos(2*pi*f*t+ (2/3)*pi - Φ),
       Theta[1] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][1],
       Theta[2] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][2],
       Theta[3] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][3],
       Theta[4] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][4],
       minusTheta[1] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][1],
       minusTheta[2] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][2],
       minusTheta[3] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][3],
       minusTheta[4] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][4],
       xdq1[1] .~ (Theta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * [x1, x2, x3])[1],
       xdq1[2] .~ (Theta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * [x1, x2, x3])[2],
       xdq2[1] .~ (minusTheta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * [x1, x2, x3])[1],
       xdq2[2] .~ (minusTheta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * [x1, x2, x3])[2],
       xdq1[1] .~ (xdq1M +  Theta * Theta * xdq2M)[1],
       xdq1[2] .~ (xdq1M +  Theta * Theta * xdq2M)[2],
       xdq2[1] .~ (xdq2M +  minusTheta * minusTheta * xdq1M)[1],
       xdq2[2] .~ (xdq2M +  minusTheta * minusTheta * xdq1M)[2],
       ]

@named ns = NonlinearSystem(eqs,[xd1M, xq1M, xd2M, xq2M, θ, Theta, minusTheta, x1, x2, x3, xdq1, xdq2], [t])
prob = NonlinearProblem(ns, ones(20), 1)

I’m using a NonlinearSystem since this is a simplified version of the problem and eventually it’ll be non linear. When running the code, I get the error ERROR: MethodError: no method matching hasmetadata(::Symbolics.Arr{Num, 2}, ::Type{Symbolics.VariableDefaultValue})

I see that this error can be fixed for ODESystem, but I don’t know what to do with a NonlinearSystem. How do I solve this?

Thanks

The issue was with the call to NonlinearSystem, I need to specify the variables one by one instead of just passing the array. I’m sure there has to be a better way to do this, but it is fixed now

@named ns = NonlinearSystem(eqs,[θ, x[1], x[2], x[3], Theta[1], Theta[2], Theta[3], Theta[4], minusTheta[1], minusTheta[2], minusTheta[3], minusTheta[4], xd1, xq1, xd2, xq2, xdq1[1], xdq1[2], xdq2[1], xdq2[2]], [t])

The code still doesn’t run, as I’m running into another error when creating the NonlinearProblem

prob = NonlinearProblem(ns, ones(20), [1])

ERROR: LoadError: MethodError: no method matching operation(::Type{Symbolics.ArrayOp{Vector{Real}}})

I just had to add collect statements so that the broadcast operator was executed. Final working code:

using Symbolics, NonlinearSolve, ModelingToolkit

V1 = 0.9
V2 = 1
V3 = 1.05
f = 50
Φ = 1
@parameters t
@variables θ Theta[1:2,1:2] minusTheta[1:2,1:2] x[1:3] xdq1[1:2] xdq2[1:2] xdq1M[1:2] xdq2M[1:2]

eqs = [
       θ ~ t,
       x[1] ~ V1 * cos(2*pi*f*t- Φ),
       x[2] ~ V2 * cos(2*pi*f*t- (2/3)*pi - Φ),
       x[3] ~ V3 * cos(2*pi*f*t+ (2/3)*pi - Φ),
       Theta[1] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][1],
       Theta[2] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][2],
       Theta[3] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][3],
       Theta[4] .~ [cos(θ) sin(θ); -sin(θ) cos(θ)][4],
       minusTheta[1] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][1],
       minusTheta[2] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][2],
       minusTheta[3] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][3],
       minusTheta[4] .~ [cos(-θ) sin(-θ); -sin(-θ) cos(-θ)][4],
       xdq1[1] .~ collect(Theta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * x)[1],
       xdq1[2] .~ collect(Theta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * x)[2],
       xdq2[1] .~ collect(minusTheta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * x)[1],
       xdq2[2] .~ collect(minusTheta * sqrt(2/3) * [1 -0.5 -0.5; 0 sqrt(3)/2 -sqrt(3)/2] * x)[2],
       xdq1[1] .~ collect(xdq1M +  Theta * Theta * xdq2M)[1],
       xdq1[2] .~ collect(xdq1M +  Theta * Theta * xdq2M)[2],
       xdq2[1] .~ collect(xdq2M +  minusTheta * minusTheta * xdq1M)[1],
       xdq2[2] .~ collect(xdq2M +  minusTheta * minusTheta * xdq1M)[2],
       ]

vars = [θ, x[1], x[2], x[3], Theta[1], Theta[2], Theta[3], Theta[4], minusTheta[1], minusTheta[2], minusTheta[3], minusTheta[4], xdq1[1], xdq1[2], xdq2[1], xdq2[2], xdq1M[1], xdq1M[2], xdq2M[1], xdq2M[2]]
@named ns = NonlinearSystem(eqs, vars, [t])
prob = NonlinearProblem(ns, ones(length(vars)), [1])
sol = solve(prob, NewtonRaphson())
1 Like

Yup, symbolic arrays still have a few issues so until that’s cleaned up collecting is a good option.

1 Like

You probably want to structural_simplify too so that way it does nonlinear tearing.