Controls Implementation using Dymola

Question

I've been having issues implementing a built-in controls algorithm (LQR) in Modelica/Dymola. I've created a simpler model that shows the problems I've been having.

Basically, I am calling LQR externally and passing it a matrix with an entry (named Ctest) that changes each time step. This Ctest is also found through an external function named findC.

The funny thing is, Dymola will run fine for if you make Ctest = 0, but if you make Ctest = 0 in an if-loop, it shows the many errors which look something like this: Unsupported: In function Modelica_LinearSystems2.Math.Matrices.dare variable AT was declared with dimension ":". That is not yet supported in dsmodel.c, and the function will fail if called in the model.

For example, I ran 3 different cases with the following code: 1. in function findC, if you define C as 0 (as copied below) everything runs fine. 2. if you instead make an if loop in which C will still be 0, the model won't simulate. I've copied the if loop below and commented it out. 3. if you keep findC as in case 1, but just uncomment "Real tether_l = 151.61;", it gives the same errors as case 2.

Any help would be greatly appreciated!

model SimplerModel
import Utilities;

Modelica.Mechanics.MultiBody.Joints.FreeMotion freeMotion(
 useQuaternions=false,
 angles_fixed=true,
r_rel_a(start={1,0,0}, fixed=true),
v_rel_a(start={0,0,0}, fixed=true),
a_rel_a(start={0,0,0}),
    angles_start={0,0,0},
    w_rel_a_fixed=true,
    w_rel_a_start={0,0,0},
    z_rel_a_fixed=false)
    annotation (Placement(transformation(extent={{-50,60},{-30,80}})));

  Modelica.Mechanics.MultiBody.Parts.BodyShape bodyShape(
    r={0,0,1},
    m=600,
    I_11=100,
    I_22=100,
    I_33=500,
    angles_start={0,0,0},
    sequence_start={1,2,3},
    w_0_start={0,0,0},
    z_0_start={0,0,0},
    r_0(start={0,0,0}),
    v_0(start={0,0,0}),
    a_0(start={0,0,0}),
    angles_fixed=false,
    w_0_fixed=false,
    z_0_fixed=false,
    r_CM={0,0,0.5})
    annotation (Placement(transformation(extent={{-10,-10},{10,10}})));
      inner Modelica.Mechanics.MultiBody.World world
    annotation (Placement(transformation(extent={{-80,60},{-60,80}})));

 Real[6,6] restMat=
 [276533.0, 0.0, 0.0, 0.0, 0.0, 0.0;
 0.0, 276533.0, 0.0, 0.0, 0.0, 0.0;
 Ctest, 0.0, 319160000.0, 0.0, 0.0, 0.0;
 0.0, 0.0, 0.0, 86086300000.0, 0.0, 0.0;
 0.0, 0.0, 0.0, 0.0, 86086300000.0, 0.0;
 0.0, 0.0, 0.0, 0.0, 0.0, 146286000.0];

 Real Ctest = Utilities.findC(bodyShape.frame_a.r_0[1]);
 Real K_cat[:,:] = Utilities.findK(restMat);

equation 
  connect(freeMotion.frame_b, bodyShape.frame_a) annotation (Line(
      points={{-30,70},{-20,70},{-20,0},{-10,0}},
      color={95,95,95},
      thickness=0.5,
      smooth=Smooth.None));
  connect(world.frame_b, freeMotion.frame_a) annotation (Line(
      points={{-60,70},{-50,70}},
      color={95,95,95},
      thickness=0.5,
      smooth=Smooth.None));
  annotation (uses(Modelica(version="3.2")), Diagram(coordinateSystem(
          preserveAspectRatio=false, extent={{-100,-100},{100,100}}), graphics));
end SimplerModel;

function findK

function findK

  import Modelica_LinearSystems2;
  input Real[6,6] restoring;

Real cyl_mass = 8.21e6;
 Real[6,6] mass = [1.253e7, 0,0,0,-2.99e8,0;
 0,1.253e7,0,2.99e8,0,0;
 0,0,1.6746e6,0,0,0;
 0,2.99e8,0,9.549e9,0,0;
 -2.99e8,0,0,0,9.549e9,0;
 0,0,0,0,0,3.4728e7];

Real[6,6] damping = [1e5,0,0,0,0,0;
0,1e5,0,0,0,0;
0,0,1.3e5,0,0,0;
0,0,0,0,0,0;
0,0,0,0,0,0;
0,0,0,0,0,1.3e7];

Real Ipitroll = 384770000;
Real Iyaw = 291440000;

protected 
Real[6,6] addMassMat = [0,0,0,0,0,0;
0,0,0,0,0,0;
0,0,cyl_mass,0,0,0;
0,0,0,Ipitroll,0,0;
0,0,0,0,Ipitroll,0;
0,0,0,0,0,Iyaw];

Real[6,6] massMat = Modelica.Math.Matrices.inv(mass + addMassMat);

Real[4, 4] A_cat = cat(1, cat(2,zeros(2,2), identity(2)), cat(2, -restoring[4:5,:]*massMat[:,4:5],-damping[4:5,:]*massMat[:,4:5]));
Real[4, 2] B_cat = cat(1, cat(1,zeros(2,2), 23/cyl_mass*identity(2)));
Real[2, 4] C_cat = cat(2, identity(2), zeros(2,2));
Real[2, 2] D_cat = zeros(2, 2);

Real[4,4] Q = [1e8,0,0,0;
 0,1e8,0,0;
 0,0,1e-8,0;
 0,0,0,1e-8];
Real[2,2] R = [1e-9,0;
 0,1e-9];

output Real K_cat[6,6];

algorithm 
K_cat := Modelica_LinearSystems2.StateSpace.Design.lqr(Modelica_LinearSystems2.StateSpace(A_cat,B_cat,C_cat,D_cat), Q, R);

end findK;

function findC

function findC

input Real x;
output Real C;

   //Real tether_l = 151.61;
   //Real slope_ForceVsHeave = 3.1928e8;
   //Real intercept_ForceVsHeave = 0;
   //Real heave = tether_l - sqrt(tether_l^2 - x^2);

algorithm 
 //if abs(x) == 0.0 then
 //C := 0;
 //else C := 0;
 //end if;
    C:=0;
end findC;

Thanks a lot! :)

Answer 1

This is the code after modifying according to the two answers (declare sizes of all matrices, protect Reals in functions, change if condition):

model SimplerModel
import OnWind.Components.Substructure.Utilities;

Modelica.Mechanics.MultiBody.Joints.FreeMotion freeMotion(
 useQuaternions=false,
 angles_fixed=true,
r_rel_a(start={0.5,0,0}, fixed=true),
v_rel_a(start={0,0,0}, fixed=true),
a_rel_a(start={0,0,0}),
    angles_start={0,0,0},
    w_rel_a_fixed=true,
    w_rel_a_start={0,0,0},
    z_rel_a_fixed=false)
    annotation (Placement(transformation(extent={{-50,60},{-30,80}})));

  Modelica.Mechanics.MultiBody.Parts.BodyShape bodyShape
    annotation (Placement(transformation(extent={{-10,-10},{10,10}})));
      inner Modelica.Mechanics.MultiBody.World world
    annotation (Placement(transformation(extent={{-80,60},{-60,80}})));

 //input Integer hi = 1;
 Real[6,6] restMat=
 [276533.0, 0.0, 0.0, 0.0, 0.0, 0.0;
 0.0, 276533.0, 0.0, 0.0, 0.0, 0.0;
 Ctest, 0.0, 319160000.0, 0.0, 0.0, 0.0;
 0.0, 0.0, 0.0, 86086300000.0, 0.0, 0.0;
 0.0, 0.0, 0.0, 0.0, 86086300000.0, 0.0;
 0.0, 0.0, 0.0, 0.0, 0.0, 146286000.0];

 Real Ctest = Utilities.findC(bodyShape.frame_b.r_0[1]);
 Real K_cat[2,4] = Utilities.findK(restMat);

equation 
  connect(freeMotion.frame_b, bodyShape.frame_a) annotation (Line(
      points={{-30,70},{-20,70},{-20,0},{-10,0}},
      color={95,95,95},
      thickness=0.5,
      smooth=Smooth.None));
  connect(world.frame_b, freeMotion.frame_a) annotation (Line(
      points={{-60,70},{-50,70}},
      color={95,95,95},
      thickness=0.5,
      smooth=Smooth.None));
  annotation (uses(Modelica(version="3.2")), Diagram(coordinateSystem(
          preserveAspectRatio=false, extent={{-100,-100},{100,100}}), graphics));
end SimplerModel;

function findK

function findK

  import Modelica_LinearSystems2;
  input Real[6,6] restoring;

protected 
Real cyl_mass = 8.21e6;
 Real[6,6] mass = [1.253e7, 0,0,0,-2.99e8,0;
 0,1.253e7,0,2.99e8,0,0;
 0,0,1.6746e6,0,0,0;
 0,2.99e8,0,9.549e9,0,0;
 -2.99e8,0,0,0,9.549e9,0;
 0,0,0,0,0,3.4728e7];

Real[6,6] damping = [1e5,0,0,0,0,0;
0,1e5,0,0,0,0;
0,0,1.3e5,0,0,0;
0,0,0,0,0,0;
0,0,0,0,0,0;
0,0,0,0,0,1.3e7];

Real Ipitroll = 384770000;
Real Iyaw = 291440000;

Real[6,6] addMassMat = [0,0,0,0,0,0;
0,0,0,0,0,0;
0,0,cyl_mass,0,0,0;
0,0,0,Ipitroll,0,0;
0,0,0,0,Ipitroll,0;
0,0,0,0,0,Iyaw];

Real[6,6] massMat = Modelica.Math.Matrices.inv(mass + addMassMat);

Real[4, 4] A_cat = cat(1, cat(2,zeros(2,2), identity(2)), cat(2, -restoring[4:5,:]*massMat[:,4:5],-damping[4:5,:]*massMat[:,4:5]));
Real[4, 2] B_cat = cat(1, cat(1,zeros(2,2), 23/cyl_mass*identity(2)));
Real[2, 4] C_cat = cat(2, identity(2), zeros(2,2));
Real[2, 2] D_cat = zeros(2, 2);

Real[4,4] Q = [1e8,0,0,0;
 0,1e8,0,0;
 0,0,1e-8,0;
 0,0,0,1e-8];
Real[2,2] R = [1e-9,0;
 0,1e-9];

output Real K_cat[2,4];

algorithm 
K_cat := Modelica_LinearSystems2.StateSpace.Design.lqr(Modelica_LinearSystems2.StateSpace(A_cat,B_cat,C_cat,D_cat), Q, R);

end findK;

function findC

function findC

input Real x;
output Real C;

   //Real tether_l = 151.61;

algorithm 
 C:=0;

/* if Modelica.Math.isEqual(abs(x),0,1e-15) then
 C := 0;
 else C := 0;
 end if;  
 */
end findC;

This should run as-is. But I've commented out necessary parts of function findC in order to make it work. In findC, uncommenting either the declaration of unused variable tether_l (which is needed in my actual code) or the if statement (also needed in my code, but written here as a trivial statement just for this simpler example) will cause it to fail.

Just to summarize, Dymola produces a long list of "Unsupported: In function __ variable __ was declared with dimension ":"". But all of these variables are actually Dymola declared in Dymola functions. Some changes cause it to work, but these changes are unrelated to the matrices.

Would appreciate any ideas or advice! :)

Thanks

Answer 2

I think Marco is right, that the problem is that the functions you are using contain variables of unknown dimensions and Dymola cannot handle this in a model. For example the dare function has:

Real AT[:, :]=transpose(A);

Try rewriting the code so that all the :'s are replaced with sizes, you may have to create functions in some cases to calculate what these sizes are going to be.

What version of Dymola are you using? I put the code, as is, into Dymola 2014 and it returned a long list of "variable ? was declared with dimension ":"."

Also I think all Real variables in functions should be protected.

Hopefully fixing this will sort out the problem.

Answer 3

Unsupported: In function Modelica_LinearSystems2.Math.Matrices.dare variable AT was declared with dimension ":"

means that you have to declare the dimension of your inputs, like:

...
input Integer n;
input Real[n] C;
...

because this, if compiled, will lead into an error:

...
input Real[:] C;
...

I think that one problem in your if-statement is that

if abs(x) == 0.0 then

should be replaced with

Modelica.Math.isEqual(abs(x),0,1e-15)

because you cannot compare the equality of two real numbers in Modelica.

I hope this helps, Marco