матрица сингулярна относительно заданной линейной системы, не разрешимой

После этого вопроса я изменил свой код на:

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;

Однако я получаю предупреждение о предварительной обработке

[1] .... Предупреждение о переводе

Итерационные переменные с нулевым начальным атрибутом по умолчанию в разорванной системе нелинейных уравнений:

     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

что не имеет смысла, но я предполагаю, что могу спокойно игнорировать и ошибку компиляции:

Матрица особенная!

недоопределенная линейная система не разрешима

о которых также ранее сообщалось здесь. если я удалю линии

Torque Tfsc;

а также

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

и изменение

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

работает отлично. Однако установка Muk в ноль, что теоретически то же самое приводит к той же ошибке, что и выше! Буду признателен, если вы поможете мне узнать, в чем проблема и как я могу ее решить.

P.S.1. В демо-версии Dymola симуляционный тест завершается без ошибок, только с предупреждением:

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. Использование JModelica, удаление noEvent и использование кода Python:

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()

он невероятно быстро решает модель для небольших значений (например, 0.1) Muk. Но снова он застревает в поисках больших ценностей. Единственные предупреждения:

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!

person Foad    schedule 05.06.2019    source источник
comment
Некоторые незначительные комментарии: это кг, а не кг, и множественное деление не допускается, поэтому кг / м / с2 следует записывать как кг / (м.с2).   -  person Hans Olsson    schedule 10.06.2019
comment
@HansOlsson Большое спасибо. Было бы здорово, если бы вы и другие ребята из Modelica были в группе Modelica Discord.   -  person Foad    schedule 13.06.2019


Ответы (2)


Вам не нужно использовать алгоритм для присвоения уравнений (даже если они находятся в цикле for и if). Я переместил их в раздел уравнений и полностью удалил ваш раздел алгоритмов:

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;

Это значительно упрощает компилятору поиск разумной сортировки и удаления сильных компонентов. Это все еще разрывается на 19, но до этого это может быть именно то, что вы ищете. Решатель Ньютона расходится после этого порога, так как я действительно не знаю, что вы здесь делаете, я, к сожалению, не могу предоставить никакого анализа результатов.

person kabdelhak    schedule 11.09.2019

Также кажется, что ваше событие, вызванное вашим if-уравнением, можно было бы чисто заменить оператором Sample. Возможно, вы захотите взглянуть на это.

person kabdelhak    schedule 11.09.2019