LinearModel.hpp

Go to the documentation of this file.

#pragma once

 
#include <icarus/core/CoreTypes.hpp>
#include <icarus/core/Error.hpp>
 
#include <Eigen/Dense>
#include <cstdio>
#include <fstream>
#include <iomanip>
#include <string>
#include <vector>
 
namespace icarus::staging {

struct LinearModel {
    Eigen::MatrixXd A; 
    Eigen::MatrixXd B; 
    Eigen::MatrixXd C; 
    Eigen::MatrixXd D; 
 
    std::vector<std::string> state_names;
    std::vector<std::string> input_names;
    std::vector<std::string> output_names;

    Eigen::VectorXd x0;
    Eigen::VectorXd u0;
    double t0 = 0.0;
 
    // =========================================================================
    // Export Methods
    // =========================================================================

    void ExportMatlab(const std::string &path) const {
        std::ofstream file(path);
        if (!file.is_open()) {
            throw IOError("Failed to open file for MATLAB export: " + path);
        }
 
        try {
            // Enable exceptions for write errors
            file.exceptions(std::ios::failbit | std::ios::badbit);
 
            file << std::setprecision(15);
            file << "% Linear state-space model\n";
            file << "% Generated by Icarus " << icarus::Version() << "\n";
            file << "% Operating point: t = " << t0 << "\n\n";
 
            // Helper to write a matrix
            auto write_matrix = [&](const std::string &name, const Eigen::MatrixXd &M) {
                file << name << " = [\n";
                for (Eigen::Index i = 0; i < M.rows(); ++i) {
                    file << "  ";
                    for (Eigen::Index j = 0; j < M.cols(); ++j) {
                        file << M(i, j);
                        if (j < M.cols() - 1)
                            file << ", ";
                    }
                    file << ";\n";
                }
                file << "];\n\n";
            };
 
            write_matrix("A", A);
            write_matrix("B", B);
            write_matrix("C", C);
            write_matrix("D", D);
 
            // Create ss object
            file << "sys = ss(A, B, C, D);\n\n";
 
            // State names
            if (!state_names.empty()) {
                file << "sys.StateName = {";
                for (std::size_t i = 0; i < state_names.size(); ++i) {
                    file << "'" << state_names[i] << "'";
                    if (i < state_names.size() - 1)
                        file << ", ";
                }
                file << "};\n";
            }
 
            // Input names
            if (!input_names.empty()) {
                file << "sys.InputName = {";
                for (std::size_t i = 0; i < input_names.size(); ++i) {
                    file << "'" << input_names[i] << "'";
                    if (i < input_names.size() - 1)
                        file << ", ";
                }
                file << "};\n";
            }
 
            // Output names
            if (!output_names.empty()) {
                file << "sys.OutputName = {";
                for (std::size_t i = 0; i < output_names.size(); ++i) {
                    file << "'" << output_names[i] << "'";
                    if (i < output_names.size() - 1)
                        file << ", ";
                }
                file << "};\n";
            }
 
            file << "\ndisp('Linear model loaded successfully.');\n";
            file.flush();
        } catch (const std::ios_base::failure &e) {
            throw IOError("Failed to write MATLAB export file '" + path + "': " + e.what());
        }
    }

    void ExportNumPy(const std::string &path) const {
        std::ofstream file(path);
        if (!file.is_open()) {
            throw IOError("Failed to open file for NumPy export: " + path);
        }
 
        try {
            // Enable exceptions for write errors
            file.exceptions(std::ios::failbit | std::ios::badbit);
 
            file << std::setprecision(15);
            file << "\"\"\"Linear state-space model\n";
            file << "Generated by Icarus " << icarus::Version() << "\n";
            file << "Operating point: t = " << t0 << "\n";
            file << "\"\"\"\n\n";
            file << "import numpy as np\n\n";
 
            // Helper to write a matrix
            auto write_matrix = [&](const std::string &name, const Eigen::MatrixXd &M) {
                file << name << " = np.array([\n";
                for (Eigen::Index i = 0; i < M.rows(); ++i) {
                    file << "    [";
                    for (Eigen::Index j = 0; j < M.cols(); ++j) {
                        file << M(i, j);
                        if (j < M.cols() - 1)
                            file << ", ";
                    }
                    file << "]";
                    if (i < M.rows() - 1)
                        file << ",";
                    file << "\n";
                }
                file << "])\n\n";
            };
 
            write_matrix("A", A);
            write_matrix("B", B);
            write_matrix("C", C);
            write_matrix("D", D);
 
            // Metadata
            file << "state_names = [";
            for (std::size_t i = 0; i < state_names.size(); ++i) {
                file << "'" << state_names[i] << "'";
                if (i < state_names.size() - 1)
                    file << ", ";
            }
            file << "]\n";
 
            file << "input_names = [";
            for (std::size_t i = 0; i < input_names.size(); ++i) {
                file << "'" << input_names[i] << "'";
                if (i < input_names.size() - 1)
                    file << ", ";
            }
            file << "]\n";
 
            file << "output_names = [";
            for (std::size_t i = 0; i < output_names.size(); ++i) {
                file << "'" << output_names[i] << "'";
                if (i < output_names.size() - 1)
                    file << ", ";
            }
            file << "]\n\n";
 
            // Operating point
            file << "x0 = np.array([";
            for (Eigen::Index i = 0; i < x0.size(); ++i) {
                file << x0(i);
                if (i < x0.size() - 1)
                    file << ", ";
            }
            file << "])\n";
 
            file << "u0 = np.array([";
            for (Eigen::Index i = 0; i < u0.size(); ++i) {
                file << u0(i);
                if (i < u0.size() - 1)
                    file << ", ";
            }
            file << "])\n";
 
            file << "t0 = " << t0 << "\n";
            file.flush();
        } catch (const std::ios_base::failure &e) {
            throw IOError("Failed to write NumPy export file '" + path + "': " + e.what());
        }
    }

    void ExportJSON(const std::string &path) const {
        std::ofstream file(path);
        if (!file.is_open()) {
            throw IOError("Failed to open file for JSON export: " + path);
        }
 
        try {
            // Enable exceptions for write errors
            file.exceptions(std::ios::failbit | std::ios::badbit);
 
            file << std::setprecision(15);
            file << "{\n";
            file << "  \"generator\": \"Icarus " << icarus::Version() << "\",\n";
            file << "  \"t0\": " << t0 << ",\n";
 
            // Helper to write a matrix
            auto write_matrix = [&](const std::string &name, const Eigen::MatrixXd &M, bool last) {
                file << "  \"" << name << "\": [\n";
                for (Eigen::Index i = 0; i < M.rows(); ++i) {
                    file << "    [";
                    for (Eigen::Index j = 0; j < M.cols(); ++j) {
                        file << M(i, j);
                        if (j < M.cols() - 1)
                            file << ", ";
                    }
                    file << "]";
                    if (i < M.rows() - 1)
                        file << ",";
                    file << "\n";
                }
                file << "  ]" << (last ? "" : ",") << "\n";
            };
 
            // Helper to write a string array
            auto write_strings = [&](const std::string &name, const std::vector<std::string> &arr,
                                     bool last) {
                file << "  \"" << name << "\": [";
                for (std::size_t i = 0; i < arr.size(); ++i) {
                    file << "\"" << arr[i] << "\"";
                    if (i < arr.size() - 1)
                        file << ", ";
                }
                file << "]" << (last ? "" : ",") << "\n";
            };
 
            // Helper to write a vector
            auto write_vector = [&](const std::string &name, const Eigen::VectorXd &v, bool last) {
                file << "  \"" << name << "\": [";
                for (Eigen::Index i = 0; i < v.size(); ++i) {
                    file << v(i);
                    if (i < v.size() - 1)
                        file << ", ";
                }
                file << "]" << (last ? "" : ",") << "\n";
            };
 
            write_strings("state_names", state_names, false);
            write_strings("input_names", input_names, false);
            write_strings("output_names", output_names, false);
            write_vector("x0", x0, false);
            write_vector("u0", u0, false);
            write_matrix("A", A, false);
            write_matrix("B", B, false);
            write_matrix("C", C, false);
            write_matrix("D", D, true);
 
            file << "}\n";
            file.flush();
        } catch (const std::ios_base::failure &e) {
            throw IOError("Failed to write JSON export file '" + path + "': " + e.what());
        }
    }
 
    // =========================================================================
    // Analysis Methods
    // =========================================================================

    [[nodiscard]] int ControllabilityRank() const {
        const Eigen::Index n = A.rows();
        const Eigen::Index m = B.cols();
 
        // If no inputs, system is not controllable
        if (m == 0 || n == 0) {
            return 0;
        }
 
        Eigen::MatrixXd W(n, n * m);
        Eigen::MatrixXd Ak = Eigen::MatrixXd::Identity(n, n);
 
        for (Eigen::Index k = 0; k < n; ++k) {
            W.block(0, k * m, n, m) = Ak * B;
            Ak = A * Ak;
        }
 
        Eigen::JacobiSVD<Eigen::MatrixXd> svd(W);
        return static_cast<int>(svd.rank());
    }

    [[nodiscard]] int ObservabilityRank() const {
        const Eigen::Index n = A.rows();
        const Eigen::Index p = C.rows();
 
        // If no outputs, system is not observable
        if (p == 0 || n == 0) {
            return 0;
        }
 
        Eigen::MatrixXd O(n * p, n);
        Eigen::MatrixXd Ak = Eigen::MatrixXd::Identity(n, n);
 
        for (Eigen::Index k = 0; k < n; ++k) {
            O.block(k * p, 0, p, n) = C * Ak;
            Ak = A * Ak;
        }
 
        Eigen::JacobiSVD<Eigen::MatrixXd> svd(O);
        return static_cast<int>(svd.rank());
    }

    [[nodiscard]] bool IsStable() const {
        Eigen::EigenSolver<Eigen::MatrixXd> solver(A, false);
        for (Eigen::Index i = 0; i < A.rows(); ++i) {
            if (solver.eigenvalues()(i).real() >= 0) {
                return false;
            }
        }
        return true;
    }

    [[nodiscard]] Eigen::VectorXcd Eigenvalues() const {
        Eigen::EigenSolver<Eigen::MatrixXd> solver(A, false);
        return solver.eigenvalues();
    }
};
 
} // namespace icarus::staging