c ++:求解非自由方程的动态系统

c++: Solving dynamic system of nonlibear equations

我正在尝试编写计算两个伺服驱动器输入转速的模型。该模型由四个方程和八个变量组成。

平均速度:v = (vd+vb)/ 2
滑滚比:srr = (vd-vb)/ v
圆盘速度:vd = 2* pi* Rd* nd
球速:vb = 2* pi* Rb* nb

在应用程序 (QT gui) 内部,模型的工作方式类似于 table,其中每个单元格代表变量。一些单元格用户填充值,其余单元格根据方程式计算。解决方案在每次输入后完成 - 如果用户插入 "vd" 和 "vb",应用程序解决 "v" 和 "srr",然后等待下一个输入。

现在,模型通过为最可能的输入组合定义的一堆 ifs 条件求解。然而,这不是一个非常优雅的解决方案,因为代码混乱且难以edit/extend。

主要problem/question是如何将这些方程式放入矩阵并求解每个输入组合(输入应该是4个变量的定义)。

//Class of the model
class Prediction : public QObject
{
    Q_OBJECT
public:
    Prediction(QObject* parent = nullptr);

    PhysicalQuantity constant(double value);

    //called to solve the model
    void evaluate();

signals:
    void callDataChanged();

private:
    const double pi = 3.141592653589793;
    int m_row;

    PhysicalQuantity m_meanSpeed  = PhysicalQuantity("v");
    PhysicalQuantity m_ballRPM = PhysicalQuantity("bRPM");
    PhysicalQuantity m_discRPM = PhysicalQuantity("dRPM");
    PhysicalQuantity m_ballSpeed = PhysicalQuantity("vb");
    PhysicalQuantity m_discSpeed = PhysicalQuantity("vd");
    PhysicalQuantity m_slideToRoll = PhysicalQuantity("srr");
    PhysicalQuantity m_discRadius = PhysicalQuantity("rd");
    PhysicalQuantity m_ballRadius = PhysicalQuantity("rb");
};
//current sollution
void Prediction::evaluate()
{
    if(this->m_meanSpeed.isEntered() && this->m_slideToRoll.isEntered() && m_discRadius.isEntered() && m_ballRadius.isEntered())
    {
        m_discSpeed.computedValue(((m_slideToRoll * m_meanSpeed + 2 * m_meanSpeed)/m_slideToRoll).value());
        m_ballSpeed.computedValue((2*m_meanSpeed - m_discSpeed).value());
        m_discRPM.computedValue((m_discSpeed / (2 * pi * m_discRadius)).value());
        m_ballRPM.computedValue((m_ballSpeed / (2 * pi * m_ballRadius)).value());
    }
...
};
//class of the variables in the model
class PhysicalQuantity
{

public:
    PhysicalQuantity();
    PhysicalQuantity(QString id);
    PhysicalQuantity(double value);
    PhysicalQuantity(const PhysicalQuantity &obj);

    //read functions
    bool isEntered() const;
    bool isConstant() const;
    double value() const;
    double matInt() const;
    QString name() const;

    //write functions
    bool setEntered(bool isEnetered);
    QString name(QString id);
    double value(double value);
    double computedValue(double value);


    //Operators
    PhysicalQuantity operator+(const PhysicalQuantity &obj);
    PhysicalQuantity operator-(const PhysicalQuantity &obj);
    PhysicalQuantity operator/(const PhysicalQuantity &obj);
    PhysicalQuantity operator=(const PhysicalQuantity &obj);
    PhysicalQuantity operator-();

protected:
    bool m_isEntered;
    bool m_isConstant;
    bool m_isComputed;
    double m_value = 1234;
    double m_defaultVal;
    QString m_id;

  • v = (vd+vb)/ 2
  • srr = (vd-vb)/ v
  • vd = 2* pi* Rd* nd
  • vb = 2* pi* Rb* nb

可能的公式是(当已经有 v = (vd+vb)/2 而不需要 srr 时,我删除了一些不需要额外变量的东西,例如 v = (vd-vb)/srr :

  • vd = 2* pi* Rd* ndvd = 2*v - vb

  • vb = 2* pi* Rb* nbvb = 2*v - vd

  • v = (vd+vb)/2

  • srr = (vd-vb)/v

  • Rd = vd/(2* pi* nd)

  • Rb = vb/(2* pi* nb)

  • nd = vd / (2* pi* Rd)

  • nb = vb/(2* pi* Rb)

所以只有10个公式,如果你喜欢按数据编程:

#include <string>
#include <vector>
#include <iostream>

const double pi = 3.141592653589793;

// index in the vector of variables
enum VARS { VD, VB, V, SSR, RD, RB, ND, NB, NVARS };

// manages a var
struct Var {
  std::string name;
  bool set;
  double value;
};

struct Formula {
  bool (*pf)(std::vector<Var> & vars); // the function to do the work
  VARS r; // the var set by the formula
  std::vector<VARS> used; // the vars used by the formula
};

// the functions doing computations, one per formula
// they are called only when the used vars are set
// return false if divide by 0

bool fvd1(std::vector<Var> & vars)
{
  vars[VD].value = 2*pi*vars[RD].value*vars[ND].value;
  return true;
}

bool fvd2(std::vector<Var> & vars)
{
  vars[VD].value = 2*vars[V].value - vars[VB].value;
  return true;
}

bool fvb1(std::vector<Var> & vars)
{
  vars[VB].value = 2*pi*vars[RB].value*vars[NB].value;
  return true;
}

bool fvb2(std::vector<Var> & vars)
{
  vars[VB].value = 2*vars[V].value - vars[VD].value;
  return true;
}

bool fv(std::vector<Var> & vars)
{
  vars[V].value = (vars[VD].value + vars[VB].value)/2;
  return true;
}

bool fssr(std::vector<Var> & vars)
{
  if (vars[V].value == 0)
    return false;

  vars[SSR].value = (vars[VD].value - vars[VB].value) / vars[V].value;
  return true;
}

bool fRd(std::vector<Var> & vars)
{
  if (vars[ND].value == 0)
    return false;

  vars[RD].value = vars[VD].value/(2 *pi * vars[ND].value);
  return true;
}

bool fRb(std::vector<Var> & vars)
{
  if (vars[NB].value == 0)
    return false;

  vars[RB].value = vars[VB].value/(2 *pi * vars[NB].value);
  return true;
}

bool fnd(std::vector<Var> & vars)
{
  if (vars[RD].value == 0)
    return false;

  vars[ND].value = vars[VD].value/(2 *pi * vars[RD].value);
  return true;
}

bool fnb(std::vector<Var> & vars)
{
  if (vars[RB].value == 0)
    return false;

  vars[NB].value = vars[VD].value/(2 *pi * vars[RB].value);
  return true;
}

const std::vector<Formula> Formulas = {
  { &fvd1, VD, { RD, ND } }, // to compute VD by fvd1 the needed vars are RD and ND
  { &fvd2, VD, { V, VB } },
  { &fvb1, VD, { RD, NB } },
  { &fvb2, VD, { V, VD } },
  { &fv, V, { VD, VB } },
  { &fssr, SSR, { VD, VB, V } },
  { &fRd, RD, { VD, ND } },
  { &fRb, RB, { VB, NB } },
  { &fnd, ND, { VD, RD } },
  { &fnb, NB, { VB, RB } }
};

// try to apply the formulas
// do in a loop in the case newly computed var(s) allows
// to set other(s) and that independently of the order 
// of the formulas
bool formulas(std::vector<Var> & vars, int & knownVars)
{
  bool modified;

  do {
    modified = false;

    for (const Formula & f : Formulas) {
      if (!vars[f.r].set) {
        bool ok = true;

        for (VARS v : f.used) {
          if (!vars[v].set) {
            ok = false;
            break;
          }
        }
        if (ok) {
          if (! (*f.pf)(vars))
            return false;
          vars[f.r].set = true;
          modified = true;
          knownVars += 1;
        }
      }
    }
  } while (modified);

  return true;
}

// to do in a terminal without graphic

int main()
{
  std::vector<Var> vars(NVARS);

  vars[VD] = { "vd", false, 0 };
  vars[VB] = { "vb", false, 0 };
  vars[V] = { "v", false, 0 };
  vars[SSR] = { "ssr", false, 0 };
  vars[RD] = { "Rd", false, 0 };
  vars[RB] = { "Rb", false, 0 };
  vars[ND] = { "nd", false, 0 };
  vars[NB] = { "nb", false, 0 };

  int knownVars = 0;

  do {
    std::cout << "enter var ( ";
    for (const Var & v : vars) {
      if (!v.set)
        std::cout << v.name << ' ';
    }
    std::cout << " ) ? ";

    std::string s;

    if (! (std::cin >> s)) {
      // EOF
      return -1;
    }

    std::vector<Var>::iterator it;

    for (it = vars.begin(); it != vars.end(); ++it) {
      if (!it->set && (it->name == s))
        break;
    }

    if (it == vars.end())
      std::cerr << "invalid name" << std::endl;
    else if (std::cout << "enter value ? ", ! (std::cin >> it->value)) {
      std::cerr << "invalid value" << std::endl;
      std::cin.clear();
      if (! (std::cin >> s)) {
        // EOF
        return -1;
      }
    }
    else {
      it->set = true;
      knownVars += 1;

      if (!formulas(vars, knownVars)) {
        std::cerr << "divide by 0, abort" << std::endl;
        return -1;
      }
    }
  } while (knownVars != NVARS);

  std::cout << "Done :";
  for (const Var & v : vars)
    std::cout << ' ' << v.name << '=' << v.value;
  std::cout << std::endl;

  return 0;
}

我不使用 Qt,而是使用标准的 C++ string/vector/IO 但是不会改变算法

编译与执行:

pi@raspberrypi:/tmp $ g++ -g -pedantic -Wextra -Wall f.cc
pi@raspberrypi:/tmp $ ./a.out
enter var ( vd vb v ssr Rd Rb nd nb  ) ? vd
enter value ? 1
enter var ( vb v ssr Rd Rb nd nb  ) ? nb
enter value ? 12
enter var ( vb v ssr Rd Rb nd  ) ? vb
enter value ? 2
enter var ( Rd nd  ) ? Rd
enter value ? 3
Done : vd=1 vb=2 v=1.5 ssr=-0.666667 Rd=3 Rb=0.0265258 nd=0.0530516 nb=12
pi@raspberrypi:/tmp $