見出し画像

リアルタイム制御向け C++行列ライブラリ


目的

行列計算をマイコンなどに実機実装したい

  • カルマンフィルタ、適応制御、モデル予測制御

しかし、世にあるライブラリは、

  • ランタイムエラーを回避できない(ゼロ割、マイナス平方根)

  • 解の精度を調整できないため、メモリ消費、計算時間が大きい

という問題があります。
そこで、リアルタイム制御に適した行列計算ライブラリを作ることにしました。
完成品は以下にあります。

Modeling-Coding-Automation-Project/python_numpy_to_cpp: リアルタイム制御向けC++行列ライブラリです。

以降、解説です。

コンセプト

  • 密行列、対角行列、スパース行列を使い分けられる

  • 反復計算の回数調整が可能

  • ゼロ割、マイナス平方根などを回避する機能を実装

以前に作成しました PythonからC++へ自動変換するツールで、Numpyの行列計算コードから簡単に置き換えられるようにインターフェースを作り込みました。

Modeling-Coding-Automation-Project/CodeLanguageConverter: 生成AIを用いて、コードの言語を自動変換します。 / Automatic conversion of the language of the code using a generative AI.

説明動画

詳細については以下の動画で説明していますのでご覧ください。

スパース行列の定義について

スパース行列はCompressed Sparse Row(CSR)方式を用いています。(別名、Compressed Row Storage、CRS形式とも言います)

一つの行列を三つの配列で表現します。

value = {1, 3, 8, 2, 4} // 値
row_index = {0, 0, 2, 1, 2} // 列インデックス
row_pointer  = {0, 1, 3, 5} // 行Headポインタ

この三つの配列で表現される元の行列は、以下のようになります。

このスパース行列を用いて、密行列と「スパース行列 * 密行列」の計算を行うと、以下のように0である要素の乗算、加算を行わないで効率的に計算できます。

Matrix<T, M, K> Y;

for (std::size_t i = 0; i < K; i++) {
  for (std::size_t j = 0; j < M; j++) {
    T sum = static_cast<T>(0);
    for (std::size_t k = A.row_pointer(j); k < A.row_pointer(j + 1); k++) {
      sum += A.value(k) * B(A.row_index(k), i);
    }
    Y(j, i) = sum;
  }
}

2024/11/04 複素数定義方法について改善

二つ目の動画の最後で、複素数の定義がbase_matrix側にあるため、使おうとすると、

#include "python_numpy.hpp"

using namespace PythonNumpy;

auto eigen_values_matrix =
      Matrix<DefDiag, Base::Matrix::Complex<double>, 3>(eigen_values_c.matrix);

のように、「Base::Matrix::Complex」と書かなければなりませんでした。上位レイヤーのPythonNumpyでは、下位レイヤーのBase::Matrixを使わないで書けるようにしたかったのですが、複素数だけはできない、というか、やろうとすると面倒でできずにいました。

いろいろ調べて分かりましたが、以下のように名前空間の置き換え機能を使えば、自然に下位レイヤーの機能を上位レイヤーでラッパーできることに気付きました。

#include "base_matrix.hpp"

namespace PythonNumpy {

template <typename T> using Complex = Base::Matrix::Complex<T>;

} // namespace PythonNumpy

これによって、以下のようにBase::Matrixを書かずに使えるようになりました。

#include "python_numpy.hpp"

using namespace PythonNumpy;

auto eigen_values_matrix =
      Matrix<DefDiag, Complex<double>, 3>(eigen_values_c.matrix);