1419 Add SECIRVVS optimal control example

cpp/CMakeLists.txt

@@ -206,6 +206,10 @@ if(MEMILIO_BUILD_BENCHMARKS)
     add_subdirectory(benchmarks)
 endif()
 
+if(MEMILIO_ENABLE_IPOPT)
+    add_subdirectory(tools/optimal_control)
+endif()
+
 if(MEMILIO_BUILD_SBML_MODELS)
     add_subdirectory(sbml_model_generation)
 endif()

cpp/tools/optimal_control/CMakeLists.txt

@@ -0,0 +1,3 @@
+if(MEMILIO_ENABLE_IPOPT)
+    add_subdirectory(secirvvs_ESID_dampings)
+endif()

cpp/tools/optimal_control/secirvvs_ESID_dampings/CMakeLists.txt

@@ -0,0 +1,33 @@
+if(MEMILIO_ENABLE_IPOPT)
+
+    # Executable for testing or running optimal control
+    add_executable(optimal_control_ESID
+        optimal_control_ESID.cpp
+        constraints/constraints.cpp
+        control_parameters/control_parameters.cpp
+        ipopt_solver/secirvvs_ipopt.cpp
+        optimization_model/optimization_model.cpp
+        optimization_settings/optimization_settings.cpp
+    )
+
+    target_link_libraries(optimal_control_ESID PRIVATE 
+        memilio 
+        ode_secirvvs 
+        Boost::boost 
+        Boost::filesystem 
+        AD::AD 
+        Ipopt::Ipopt
+    )
+
+    target_include_directories(optimal_control_ESID PRIVATE
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/constraints>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/control_parameters>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/helpers>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/ipopt_solver>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/optimization_model>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/optimization_settings>
+        $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+    )
+
+endif()

cpp/tools/optimal_control/secirvvs_ESID_dampings/check_feasibility.h

@@ -0,0 +1,119 @@
+#pragma once
+
+#include "optimization_settings/optimization_settings.h"
+#include "control_parameters/damping_controls.h"
+#include "constraints/update_constraints.h"
+#include "helpers/integrator_selector.h"
+#include "helpers/make_time_grid.h"
+
+#include <cstddef>
+#include <filesystem>
+#include <memory>
+#include <vector>
+#include <fstream>
+#include <stdexcept>
+#include <iomanip>
+
+void check_constraint_feasibility(const SecirvvsOptimization& settings)
+{
+    // Set most restrictive controls to validate if the constraints allow for a feasible solution.
+    std::vector<double> parameters(settings.num_control_parameters() * settings.num_control_intervals());
+    for (size_t interval = 0; interval < settings.num_control_intervals(); ++interval) {
+        for (size_t i = 0; i < settings.num_control_parameters(); ++i) {
+            size_t index      = interval * settings.num_control_parameters() + i;
+            parameters[index] = settings.control_parameters()[i].max();
+        }
+    }
+
+    // Simulate with these controls and evaluate the constraints.
+    std::vector<double> path_constraint_values(settings.num_path_constraints(), 0.0);
+    std::vector<double> terminal_constraint_values(settings.num_terminal_constraints(), 0.0);
+
+    auto graph_model = settings.optimization_model().get_graph_model<double>();
+
+    // Set the control parameters for each node simulation
+    for (auto& node : graph_model.nodes()) {
+        set_control_dampings<double>(settings, node.property.get_simulation().get_model(), parameters);
+    }
+    // Set the integrator for each node simulation
+    for (auto& node : graph_model.nodes()) {
+        node.property.get_simulation().set_integrator_core(
+            std::move(make_integrator<double>(settings.integrator_type(), settings.dt())));
+    }
+
+    // Create graph simulation
+    auto graph_simulation = mio::make_mobility_sim<double>(settings.t0(), settings.dt(), std::move(graph_model));
+
+    // Step through the time grid and evaluate the path constraints
+    std::vector<double> time_steps = make_time_grid<double>(settings.t0(), settings.tmax(), settings.num_intervals());
+
+    // Simulate and update constraints
+    update_path_constraint<double>(settings, graph_simulation.get_graph(), path_constraint_values);
+    for (size_t interval = 0; interval < settings.num_intervals(); interval++) {
+        graph_simulation.advance(time_steps[interval + 1]);
+        update_path_constraint<double>(settings, graph_simulation.get_graph(), path_constraint_values);
+    }
+    update_terminal_constraint<double>(settings, graph_simulation.get_graph(), terminal_constraint_values);
+
+    // Print out the path constraint values
+    std::cout << "Lower bound for path constraints:\n";
+    for (size_t i = 0; i < settings.num_path_constraints(); ++i) {
+        const Constraint& constraint = settings.path_constraints()[i];
+        std::cout << "  [" << i << "] \"" << constraint.name() << "\": " << path_constraint_values[i] << "\n";
+    }
+    // Print out the terminal constraint values
+    std::cout << "Lower bound on terminal constraints:\n";
+    for (size_t i = 0; i < settings.num_terminal_constraints(); ++i) {
+        const Constraint& constraint = settings.terminal_constraints()[i];
+        std::cout << "  [" << i << "] \"" << constraint.name() << "\": " << terminal_constraint_values[i] << "\n";
+    }
+
+    // ---------------------------- //
+    // Check constraint feasibility //
+    bool violated = false;
+    std::ostringstream report;
+
+    report << std::fixed << std::setprecision(6);
+    report << "Constraint feasibility report\n";
+    report << "=====================================\n\n";
+
+    // --- Path constraints (check only upper bounds) ---
+    for (size_t i = 0; i < settings.num_path_constraints(); ++i) {
+        const auto& constraint = settings.path_constraints()[i];
+        double value           = path_constraint_values[i];
+        double max_v           = constraint.max();
+
+        if (value > max_v) {
+            violated = true;
+            report << "Path constraint violated: [" << i << "] \"" << constraint.name() << "\" value = " << value
+                   << " exceeds upper bound " << max_v << "\n";
+        }
+    }
+
+    // --- Terminal constraints (check only upper bounds) ---
+    for (size_t i = 0; i < settings.num_terminal_constraints(); ++i) {
+        const auto& constraint = settings.terminal_constraints()[i];
+        double value           = terminal_constraint_values[i];
+        double max_v           = constraint.max();
+
+        if (value > max_v) {
+            violated = true;
+            report << "Terminal constraint violated: [" << i << "] \"" << constraint.name() << "\" value = " << value
+                   << " exceeds upper bound " << max_v << "\n";
+        }
+    }
+
+    // ---------------------------------- //
+    // Report and throw error if violated //
+    if (violated) {
+        std::ofstream log_file("constraint_violation.log");
+        log_file << report.str();
+        log_file.close();
+
+        std::cerr << "Constraint feasibility check failed. See 'constraint_violation.log' for details.\n";
+        throw std::runtime_error("Constraint feasibility check failed. See 'constraint_violation.log'.");
+    }
+    else {
+        std::cout << "All constraints are feasible.\n";
+    }
+}

cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.cpp

@@ -0,0 +1,48 @@
+#include "constraints.h"
+
+Constraint::Constraint(std::string name, std::pair<double, double> allowed_range, std::optional<size_t> node_id)
+    : m_name(std::move(name))
+    , m_range(allowed_range)
+    , m_node_id(node_id)
+{
+}
+
+const std::string& Constraint::name() const
+{
+    return m_name;
+}
+
+std::pair<double, double> Constraint::range() const
+{
+    return m_range;
+}
+
+double Constraint::min() const
+{
+    return m_range.first;
+}
+
+double Constraint::max() const
+{
+    return m_range.second;
+}
+
+std::optional<size_t> Constraint::node_id() const
+{
+    return m_node_id;
+}
+
+void Constraint::set_name(const std::string& new_name)
+{
+    m_name = new_name;
+}
+
+void Constraint::set_range(const std::pair<double, double>& new_range)
+{
+    m_range = new_range;
+}
+
+void Constraint::set_node_id(std::optional<size_t> new_node_id)
+{
+    m_node_id = new_node_id;
+}

cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.h

@@ -0,0 +1,26 @@
+#pragma once
+
+#include <string>
+#include <utility>
+#include <optional>
+
+class Constraint
+{
+public:
+    Constraint(std::string name, std::pair<double, double> allowed_range, std::optional<size_t> node_id);
+
+    const std::string& name() const;
+    std::pair<double, double> range() const;
+    double min() const;
+    double max() const;
+    std::optional<size_t> node_id() const;
+
+    void set_name(const std::string& new_name);
+    void set_range(const std::pair<double, double>& new_range);
+    void set_node_id(std::optional<size_t> new_node_id);
+
+private:
+    std::string m_name;
+    std::pair<double, double> m_range;
+    std::optional<size_t> m_node_id;
+};