2

Following this question, I modified my code to:

model test
  // types
  type Mass = Real(unit = "Kg", min = 0);
  type Length = Real(unit = "m");
  type Area = Real(unit = "m2", min = 0);
  type Force = Real(unit = "Kg.m/s2");
  type Pressure = Real(unit = "Kg/m/s2");
  type Torque = Real(unit = "Kg.m2/s2");
  type Velocity = Real(unit = "m/s");
  type Time = Real(unit = "s");

  // constants
  constant Real pi = 2 * Modelica.Math.asin(1.0);
  parameter Mass Mp = 0.01;
  parameter Length r1 = 0.010;
  parameter Length r3 = 0.004;
  parameter Integer n = 3;
  parameter Area A = 0.020 * 0.015;
  parameter Time Stepping = 1.0;
  parameter Real DutyCycle = 1.0;
  parameter Pressure Pin = 500000;
  parameter Real Js = 1;
  //parameter Real Muk = 0.0;
  parameter Real Muk = 0.158;

  // variables
  Length x[n];
  Velocity vx[n];
  Real theta;
  Real vt;
  Pressure P[n];
  Force Fnsp[n];
  Torque Tfsc;

initial equation
  theta = 0;
  vt = 0;

algorithm
  for i in 1:n loop
    if noEvent((i - 1) * Stepping < mod(time, n * Stepping)) and noEvent(mod(time, n * Stepping) < Stepping * ((i - 1) + DutyCycle)) then
      P[i] := Pin;
    else
      P[i] := 0;
    end if;
  end for;
  Tfsc := -r3 * Muk * sign(vt) * abs(sum(Fnsp));

equation
  vx = der(x);
  vt = der(theta);
  x = r1 * {sin(theta + (i - 1) * 2 * pi / n) for i in 1:n};
  Mp * der(vx) + P * A = Fnsp;
  Js * der(theta) = Tfsc - r1 * Fnsp * {cos(theta + (i - 1) * 2 * pi / n) for i in 1:n};
  // Js * der(theta) = - r1 * Fnsp * {cos(theta + (i - 1) * 2 * pi / n) for i in 1:n};

  annotation(
    experiment(StartTime = 0, StopTime = 30, Tolerance = 1e-06, Interval = 0.03),
    __OpenModelica_simulationFlags(lv = "LOG_STATS", outputFormat = "mat", s = "dassl"));
end test;

However, I get the preprocessing warning of

[1] .... Translation Warning

Iteration variables with default zero start attribute in torn nonlinear equation system:

     Fnsp[3]:VARIABLE(unit = "Kg.m/s2" )  type: Real  [3]
     Fnsp[2]:VARIABLE(unit = "Kg.m/s2" )  type: Real  [3]
     Fnsp[1]:VARIABLE(unit = "Kg.m/s2" )  type: Real  [3]
     $DER.vt:VARIABLE()  type: Real

which doesn't make sense but I assume I can safely ignore, and the compiling error of:

Matrix singular!

under-determined linear system not solvable

which had also been previously reported here. if I remove the lines

Torque Tfsc;

and

Tfsc := -r3 * Muk * sign(vt) * abs(sum(Fnsp));

and changing

Js * der(theta) = - r1 * Fnsp * {cos(theta + (i - 1) * 2 * pi / n) for i in 1:n};

works perfectly fine. However, setting Muk to zero, which theoretically the same thing leads to the same error as above! I would appreciate if you could help me know what is the problem and how I can resolve it.

P.S.1. On the demo version of Dymola the simulation test finishes with no errors, only the warning:

Some variables are iteration variables of the initialization problem:
but they are not given any explicit start values. Zero will be used.
Iteration variables:
der(theta, 2)
P[1]
P[2]
P[3]

P.S.2. Using JModelica, removing the noEvent and using the python code:

model_name = 'test'
mo_file = 'test.mo'

from pymodelica import compile_fmu
from pyfmi import load_fmu

my_fmu = compile_fmu(model_name, mo_file)

myModel = load_fmu('test.fmu')
res = myModel.simulate(final_time=30)

theta = res['theta']
t = res['time']

import matplotlib.pyplot as plt
plt.plot(t, theta)
plt.show()

it solves the model blazingly fast for small values (e.g. 0.1) of Muk. But again it gets stuck for bigger values. The only warnings are:

Warning at line 30, column 3, in file 'test.mo':
  Iteration variable "Fnsp[2]" is missing start value!
Warning at line 30, column 3, in file 'test.mo':
  Iteration variable "Fnsp[3]" is missing start value!
Warning in flattened model:
  Iteration variable "der(_der_theta)" is missing start value!
Foad S. Farimani
  • 12,396
  • 15
  • 78
  • 193
  • 1
    Some minor comments: it's "kg" not "Kg" and multiple divisions are not allowed so "Kg/m/s2" should be written as "kg/(m.s2)". – Hans Olsson Jun 10 '19 at 12:38
  • @HansOlsson Thanks a lot. It would be great to have you and other Modelica guys on [Modelica Discord group](https://discordapp.com/invite/bp2yeYU). – Foad S. Farimani Jun 13 '19 at 06:34

2 Answers2

1

You do not need to use algorithm for the assignments of equations (even if they are in a for-loop and an if). I moved them to the equation section and removed your algorithm section completely:

model test
  // types
  type Mass = Real(unit = "Kg", min = 0);
  type Length = Real(unit = "m");
  type Area = Real(unit = "m2", min = 0);
  type Force = Real(unit = "Kg.m/s2");
  type Pressure = Real(unit = "Kg/m/s2");
  type Torque = Real(unit = "Kg.m2/s2");
  type Velocity = Real(unit = "m/s");
  type Time = Real(unit = "s");

  // constants
  constant Real pi = 2 * Modelica.Math.asin(1.0);
  parameter Mass Mp = 0.01;
  parameter Length r1 = 0.010;
  parameter Length r3 = 0.004;
  parameter Integer n = 3;
  parameter Area A = 0.020 * 0.015;
  parameter Time Stepping = 1.0;
  parameter Real DutyCycle = 1.0;
  parameter Pressure Pin = 500000;
  parameter Real Js = 1;
  //parameter Real Muk = 0.0;
  parameter Real Muk = 0.158;

  // variables
  Length x[n];
  Velocity vx[n];
  Real theta;
  Real vt;
  Pressure P[n];
  Force Fnsp[n];
  Torque Tfsc;

initial equation
  theta = 0;
  vt = 0;
equation

  for i in 1:n loop
    if noEvent((i - 1) * Stepping < mod(time, n * Stepping)) and noEvent(mod(time, n * Stepping) < Stepping * ((i - 1) + DutyCycle)) then
      P[i] = Pin;
    else
      P[i] = 0;
    end if;
  end for;
  Tfsc = -r3 * Muk * sign(vt) * abs(sum(Fnsp));

  vx = der(x);
  vt = der(theta);
  x = r1 * {sin(theta + (i - 1) * 2 * pi / n) for i in 1:n};
  Mp * der(vx) + P * A = Fnsp;
  Js * der(theta) = Tfsc - r1 * Fnsp * {cos(theta + (i - 1) * 2 * pi / n) for i in 1:n};
  // Js * der(theta) = - r1 * Fnsp * {cos(theta + (i - 1) * 2 * pi / n) for i in 1:n};
end test;

This makes it far easier for the compiler to find a sensible sorting and tearing for the strong components. This still breaks at 19s but before that it may just be what you are looking for. The newton solver diverges after that threshold, since i don't really know what you are doing here i unfortunately cannot provide any analysis of the results.

kabdelhak
  • 687
  • 3
  • 8
1

Also it seems like your event triggered by your if-equation could be cleanly replaced by a Sample operator. You might want to have a look at that.

kabdelhak
  • 687
  • 3
  • 8