diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 3ce345b3fd..6bb0763848 100644 --- a/cpp/CMakeLists.txt +++ b/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() diff --git a/cpp/tools/optimal_control/CMakeLists.txt b/cpp/tools/optimal_control/CMakeLists.txt new file mode 100644 index 0000000000..1e79fac971 --- /dev/null +++ b/cpp/tools/optimal_control/CMakeLists.txt @@ -0,0 +1,3 @@ +if(MEMILIO_ENABLE_IPOPT) + add_subdirectory(secirvvs_ESID_dampings) +endif() diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/CMakeLists.txt b/cpp/tools/optimal_control/secirvvs_ESID_dampings/CMakeLists.txt new file mode 100644 index 0000000000..39f552e01d --- /dev/null +++ b/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 + $ + $ + $ + $ + $ + $ + $ + $ + ) + +endif() diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/check_feasibility.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/check_feasibility.h new file mode 100644 index 0000000000..019bdb3024 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include + +void check_constraint_feasibility(const SecirvvsOptimization& settings) +{ + // Set most restrictive controls to validate if the constraints allow for a feasible solution. + std::vector 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 path_constraint_values(settings.num_path_constraints(), 0.0); + std::vector terminal_constraint_values(settings.num_terminal_constraints(), 0.0); + + auto graph_model = settings.optimization_model().get_graph_model(); + + // Set the control parameters for each node simulation + for (auto& node : graph_model.nodes()) { + set_control_dampings(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(settings.integrator_type(), settings.dt()))); + } + + // Create graph simulation + auto graph_simulation = mio::make_mobility_sim(settings.t0(), settings.dt(), std::move(graph_model)); + + // Step through the time grid and evaluate the path constraints + std::vector time_steps = make_time_grid(settings.t0(), settings.tmax(), settings.num_intervals()); + + // Simulate and update constraints + update_path_constraint(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(settings, graph_simulation.get_graph(), path_constraint_values); + } + update_terminal_constraint(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"; + } +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.cpp b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.cpp new file mode 100644 index 0000000000..ac1d55f54b --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.cpp @@ -0,0 +1,48 @@ +#include "constraints.h" + +Constraint::Constraint(std::string name, std::pair allowed_range, std::optional 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 Constraint::range() const +{ + return m_range; +} + +double Constraint::min() const +{ + return m_range.first; +} + +double Constraint::max() const +{ + return m_range.second; +} + +std::optional 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& new_range) +{ + m_range = new_range; +} + +void Constraint::set_node_id(std::optional new_node_id) +{ + m_node_id = new_node_id; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.h new file mode 100644 index 0000000000..f559d98e80 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/constraints.h @@ -0,0 +1,26 @@ +#pragma once + +#include +#include +#include + +class Constraint +{ +public: + Constraint(std::string name, std::pair allowed_range, std::optional node_id); + + const std::string& name() const; + std::pair range() const; + double min() const; + double max() const; + std::optional node_id() const; + + void set_name(const std::string& new_name); + void set_range(const std::pair& new_range); + void set_node_id(std::optional new_node_id); + +private: + std::string m_name; + std::pair m_range; + std::optional m_node_id; +}; diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/infection_state_utils.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/infection_state_utils.h new file mode 100644 index 0000000000..60bf8d31a8 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/infection_state_utils.h @@ -0,0 +1,167 @@ +#pragma once + +#include "models/ode_secirvvs/model.h" + +#include +#include +#include +#include +#include +#include + +inline size_t num_infection_states() +{ + return static_cast(mio::osecirvvs::InfectionState::Count); +} + +inline std::string infection_state_to_string(mio::osecirvvs::InfectionState state) +{ + using InfectionState = mio::osecirvvs::InfectionState; + + switch (state) { + case InfectionState::SusceptibleNaive: + return "SusceptibleNaive"; + case InfectionState::SusceptiblePartialImmunity: + return "SusceptiblePartialImmunity"; + case InfectionState::SusceptibleImprovedImmunity: + return "SusceptibleImprovedImmunity"; + case InfectionState::ExposedNaive: + return "ExposedNaive"; + case InfectionState::ExposedPartialImmunity: + return "ExposedPartialImmunity"; + case InfectionState::ExposedImprovedImmunity: + return "ExposedImprovedImmunity"; + case InfectionState::InfectedNoSymptomsNaive: + return "InfectedNoSymptomsNaive"; + case InfectionState::InfectedNoSymptomsPartialImmunity: + return "InfectedNoSymptomsPartialImmunity"; + case InfectionState::InfectedNoSymptomsImprovedImmunity: + return "InfectedNoSymptomsImprovedImmunity"; + case InfectionState::InfectedNoSymptomsNaiveConfirmed: + return "InfectedNoSymptomsNaiveConfirmed"; + case InfectionState::InfectedNoSymptomsPartialImmunityConfirmed: + return "InfectedNoSymptomsPartialImmunityConfirmed"; + case InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed: + return "InfectedNoSymptomsImprovedImmunityConfirmed"; + case InfectionState::InfectedSymptomsNaive: + return "InfectedSymptomsNaive"; + case InfectionState::InfectedSymptomsPartialImmunity: + return "InfectedSymptomsPartialImmunity"; + case InfectionState::InfectedSymptomsImprovedImmunity: + return "InfectedSymptomsImprovedImmunity"; + case InfectionState::InfectedSymptomsNaiveConfirmed: + return "InfectedSymptomsNaiveConfirmed"; + case InfectionState::InfectedSymptomsPartialImmunityConfirmed: + return "InfectedSymptomsPartialImmunityConfirmed"; + case InfectionState::InfectedSymptomsImprovedImmunityConfirmed: + return "InfectedSymptomsImprovedImmunityConfirmed"; + case InfectionState::InfectedSevereNaive: + return "InfectedSevereNaive"; + case InfectionState::InfectedSeverePartialImmunity: + return "InfectedSeverePartialImmunity"; + case InfectionState::InfectedSevereImprovedImmunity: + return "InfectedSevereImprovedImmunity"; + case InfectionState::InfectedCriticalNaive: + return "InfectedCriticalNaive"; + case InfectionState::InfectedCriticalPartialImmunity: + return "InfectedCriticalPartialImmunity"; + case InfectionState::InfectedCriticalImprovedImmunity: + return "InfectedCriticalImprovedImmunity"; + case InfectionState::DeadNaive: + return "DeadNaive"; + case InfectionState::DeadPartialImmunity: + return "DeadPartialImmunity"; + case InfectionState::DeadImprovedImmunity: + return "DeadImprovedImmunity"; + case InfectionState::Count: + return "Count"; + default: + return "InfectionState Unknown"; + } +} + +inline std::vector query_infection_states(const std::string& state_query) +{ + using InfectionState = mio::osecirvvs::InfectionState; + + static const std::vector> all_states = { + {"SusceptibleNaive", InfectionState::SusceptibleNaive}, + {"SusceptiblePartialImmunity", InfectionState::SusceptiblePartialImmunity}, + {"SusceptibleImprovedImmunity", InfectionState::SusceptibleImprovedImmunity}, + {"ExposedNaive", InfectionState::ExposedNaive}, + {"ExposedPartialImmunity", InfectionState::ExposedPartialImmunity}, + {"ExposedImprovedImmunity", InfectionState::ExposedImprovedImmunity}, + {"InfectedNoSymptomsNaive", InfectionState::InfectedNoSymptomsNaive}, + {"InfectedNoSymptomsPartialImmunity", InfectionState::InfectedNoSymptomsPartialImmunity}, + {"InfectedNoSymptomsImprovedImmunity", InfectionState::InfectedNoSymptomsImprovedImmunity}, + {"InfectedNoSymptomsNaiveConfirmed", InfectionState::InfectedNoSymptomsNaiveConfirmed}, + {"InfectedNoSymptomsPartialImmunityConfirmed", InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}, + {"InfectedNoSymptomsImprovedImmunityConfirmed", InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}, + {"InfectedSymptomsNaive", InfectionState::InfectedSymptomsNaive}, + {"InfectedSymptomsPartialImmunity", InfectionState::InfectedSymptomsPartialImmunity}, + {"InfectedSymptomsImprovedImmunity", InfectionState::InfectedSymptomsImprovedImmunity}, + {"InfectedSymptomsNaiveConfirmed", InfectionState::InfectedSymptomsNaiveConfirmed}, + {"InfectedSymptomsPartialImmunityConfirmed", InfectionState::InfectedSymptomsPartialImmunityConfirmed}, + {"InfectedSymptomsImprovedImmunityConfirmed", InfectionState::InfectedSymptomsImprovedImmunityConfirmed}, + {"InfectedSevereNaive", InfectionState::InfectedSevereNaive}, + {"InfectedSeverePartialImmunity", InfectionState::InfectedSeverePartialImmunity}, + {"InfectedSevereImprovedImmunity", InfectionState::InfectedSevereImprovedImmunity}, + {"InfectedCriticalNaive", InfectionState::InfectedCriticalNaive}, + {"InfectedCriticalPartialImmunity", InfectionState::InfectedCriticalPartialImmunity}, + {"InfectedCriticalImprovedImmunity", InfectionState::InfectedCriticalImprovedImmunity}, + {"DeadNaive", InfectionState::DeadNaive}, + {"DeadPartialImmunity", InfectionState::DeadPartialImmunity}, + {"DeadImprovedImmunity", InfectionState::DeadImprovedImmunity}}; + + // Split query into separate subqueries + std::vector and_group_tokens; + std::vector or_group_tokens; + + std::istringstream iss(state_query); + std::string word; + while (iss >> word) { + if (word.find('+') != std::string::npos) { + std::istringstream substream(word); + std::string sub; + while (std::getline(substream, sub, '+')) { + and_group_tokens.push_back(sub); + } + } + else { + or_group_tokens.push_back(word); + } + } + + std::vector result; + + if (!and_group_tokens.empty()) { + for (const auto& [n, state] : all_states) { + const auto& name = n; // create a named reference or copy + bool match_all = + std::all_of(and_group_tokens.begin(), and_group_tokens.end(), [&](const std::string& token) { + return name.find(token) != std::string::npos; + }); + if (match_all) { + result.push_back(state); + } + } + } + + if (!or_group_tokens.empty()) { + for (const auto& [n, state] : all_states) { + const auto& name = n; + bool match_any = std::any_of(or_group_tokens.begin(), or_group_tokens.end(), [&](const std::string& token) { + return name.find(token) != std::string::npos; + }); + if (match_any && std::find(result.begin(), result.end(), state) == result.end()) { + result.push_back(state); + } + } + } + + if (result.empty()) { + throw std::invalid_argument("No matching InfectionState for query: " + state_query); + } + + return result; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/update_constraints.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/update_constraints.h new file mode 100644 index 0000000000..d6b6899846 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/constraints/update_constraints.h @@ -0,0 +1,85 @@ +#pragma once + +#include "constraints.h" +#include "infection_state_utils.h" + +#include "../optimization_settings/optimization_settings.h" + +#include +#include +#include +#include + +template +void update_path_constraint(const SecirvvsOptimization& settings, + mio::Graph, mio::MobilityEdge>& graph, + std::vector& path_constraint_values) +{ + assert(path_constraint_values.size() == settings.num_path_constraints()); + using std::max; + + for (size_t i = 0; i < settings.num_path_constraints(); i++) { + const Constraint& constraint = settings.path_constraints()[i]; + auto states = query_infection_states(constraint.name()); + std::optional node_id = constraint.node_id(); + + FP value = 0.0; + for (size_t node_idx = 0; node_idx < graph.nodes().size(); node_idx++) { + const auto& node = graph.nodes()[node_idx].property; + const auto& node_simulation = node.get_simulation(); + const auto& node_model = node_simulation.get_model(); + size_t graph_node_id = graph.nodes()[node_idx].id; + size_t num_age_groups = static_cast(node_model.parameters.get_num_groups()); + + const mio::TimeSeries& simulation_result = node_simulation.get_result(); + Eigen::Ref> latest_state = simulation_result.get_last_value(); + + if (!node_id.has_value() || node_id.value() == graph_node_id) { + for (size_t age_group = 0; age_group < num_age_groups; age_group++) { + for (auto state : states) { + size_t index = age_group * num_infection_states() + static_cast(state); + value += latest_state[index]; + } + } + } + } + path_constraint_values[i] = max(path_constraint_values[i], value); + } +} + +template +void update_terminal_constraint(const SecirvvsOptimization& settings, + mio::Graph, mio::MobilityEdge>& graph, + std::vector& terminal_constraint_values) +{ + assert(terminal_constraint_values.size() == settings.num_terminal_constraints()); + using std::max; + + for (size_t i = 0; i < settings.num_terminal_constraints(); i++) { + const Constraint& constraint = settings.terminal_constraints()[i]; + auto states = query_infection_states(constraint.name()); + std::optional node_id = constraint.node_id(); + + FP value = 0.0; + for (size_t node_idx = 0; node_idx < graph.nodes().size(); node_idx++) { + const auto& node = graph.nodes()[node_idx].property; + const auto& node_simulation = node.get_simulation(); + const auto& node_model = node_simulation.get_model(); + size_t graph_node_id = graph.nodes()[node_idx].id; + size_t num_age_groups = static_cast(node_model.parameters.get_num_groups()); + + const mio::TimeSeries& simulation_result = node_simulation.get_result(); + Eigen::Ref> latest_state = simulation_result.get_last_value(); + + if (!node_id.has_value() || node_id.value() == graph_node_id) { + for (size_t age_group = 0; age_group < num_age_groups; age_group++) { + for (auto state : states) { + size_t index = age_group * num_infection_states() + static_cast(state); + value += latest_state[index]; + } + } + } + } + terminal_constraint_values[i] = value; + } +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_activation.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_activation.h new file mode 100644 index 0000000000..10831a963b --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_activation.h @@ -0,0 +1,47 @@ +#pragma once + +#include + +enum class ControlActivation +{ + Linear, + Sigmoid +}; + +struct ControlActivationFunction { + ControlActivationFunction(ControlActivation activation) + : m_activation(activation) + { + } + + template + FP operator()(FP x) const + { + switch (m_activation) { + case ControlActivation::Linear: + return linear(x); + case ControlActivation::Sigmoid: + return sigmoid(x); + } + return x; + } + +private: + ControlActivation m_activation; + + // Linear activation: f(x) = x + template + static FP linear(FP x) + { + return x; + } + + // Sigmoid activation: f(x) = 0.5 + 0.5 * tanh(gain*(x-0.5)) / tanh(0.5 * gain) + template + static FP sigmoid(FP x) + { + using std::tanh; + FP gain = 100.0; + return 0.5 + 0.5 * tanh(gain * (x - 0.5)) / tanh(0.5 * gain); + } +}; diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_parameters.cpp b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_parameters.cpp new file mode 100644 index 0000000000..ce96a4e4bd --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_parameters.cpp @@ -0,0 +1,118 @@ +#include "control_parameters.h" + +DampingControlParameter::DampingControlParameter(const std::string& name, std::pair allowed_range, + double effectiveness_value, double cost_value, mio::DampingLevel level, + mio::DampingType type, std::vector locations, + const Eigen::VectorXd& group_weights) + : m_name(name) + , m_range(std::move(allowed_range)) + , m_effectiveness(effectiveness_value) + , m_cost(cost_value) + , m_level(std::move(level)) + , m_type(std::move(type)) + , m_locations(std::move(locations)) + , m_group_weights(group_weights) +{ + if (m_range.first > m_range.second) { + throw std::invalid_argument("Invalid range: min must not exceed max."); + } + if (m_group_weights.size() == 0) { + throw std::invalid_argument("Group weights vector must not be empty."); + } +} + +std::string DampingControlParameter::name() const +{ + return m_name; +} + +std::pair DampingControlParameter::range() const +{ + return m_range; +} + +double DampingControlParameter::min() const +{ + return m_range.first; +} + +double DampingControlParameter::max() const +{ + return m_range.second; +} + +double DampingControlParameter::effectiveness() const +{ + return m_effectiveness; +} + +double DampingControlParameter::cost() const +{ + return m_cost; +} + +const mio::DampingLevel& DampingControlParameter::level() const +{ + return m_level; +} + +const mio::DampingType& DampingControlParameter::type() const +{ + return m_type; +} + +const std::vector& DampingControlParameter::locations() const +{ + return m_locations; +} + +const Eigen::VectorXd& DampingControlParameter::group_weights() const +{ + return m_group_weights; +} + +void DampingControlParameter::set_name(const std::string& new_name) +{ + m_name = new_name; +} + +void DampingControlParameter::set_range(const std::pair& new_range) +{ + if (new_range.first > new_range.second) { + throw std::invalid_argument("Invalid range: min must not exceed max."); + } + m_range = new_range; +} + +void DampingControlParameter::set_effectiveness(double new_effectiveness) +{ + m_effectiveness = new_effectiveness; +} + +void DampingControlParameter::set_cost(double new_cost) +{ + m_cost = new_cost; +} + +void DampingControlParameter::set_level(const mio::DampingLevel& new_level) +{ + m_level = new_level; +} + +void DampingControlParameter::set_type(const mio::DampingType& new_type) +{ + m_type = new_type; +} + +void DampingControlParameter::set_locations(const std::vector& new_locations) +{ + m_locations = new_locations; +} + +void DampingControlParameter::set_group_weights(const Eigen::VectorXd& new_group_weights) +{ + if (new_group_weights.size() == 0) { + throw std::invalid_argument("Group weights vector must not be empty."); + } + m_group_weights = new_group_weights; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_parameters.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_parameters.h new file mode 100644 index 0000000000..3914c21099 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/control_parameters.h @@ -0,0 +1,59 @@ +#pragma once + +#include "memilio/math/eigen.h" +#include "memilio/epidemiology/damping.h" +#include "memilio/epidemiology/damping_sampling.h" + +#include +#include +#include + +class DampingControlParameter +{ +public: + DampingControlParameter(const std::string& name, std::pair allowed_range, + double effectiveness_value, double cost_value, mio::DampingLevel level, + mio::DampingType type, std::vector locations, const Eigen::VectorXd& group_weights); + + std::string name() const; + std::pair range() const; + double min() const; + double max() const; + double effectiveness() const; + double cost() const; + const mio::DampingLevel& level() const; + const mio::DampingType& type() const; + const std::vector& locations() const; + const Eigen::VectorX& group_weights() const; + + void set_name(const std::string& new_name); + void set_range(const std::pair& new_range); + void set_effectiveness(double new_effectiveness); + void set_cost(double new_cost); + void set_level(const mio::DampingLevel& new_level); + void set_type(const mio::DampingType& new_type); + void set_locations(const std::vector& new_locations); + void set_group_weights(const Eigen::VectorX& new_group_weights); + + template + mio::DampingSampling make_damping_sampling(const mio::UncertainValue& strength, + mio::SimulationTime time) const; + +private: + std::string m_name; + std::pair m_range; + double m_effectiveness; + double m_cost; + mio::DampingLevel m_level; + mio::DampingType m_type; + std::vector m_locations; + Eigen::VectorX m_group_weights; +}; + +template +mio::DampingSampling DampingControlParameter::make_damping_sampling(const mio::UncertainValue& strength, + mio::SimulationTime time) const +{ + return mio::DampingSampling(m_effectiveness * strength, m_level, m_type, time, m_locations, + m_group_weights.cast()); +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/damping_controls.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/damping_controls.h new file mode 100644 index 0000000000..79915bedc3 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/control_parameters/damping_controls.h @@ -0,0 +1,39 @@ +#pragma once + +#include "control_parameters.h" + +#include "models/ode_secirvvs/model.h" + +#include "../optimization_settings/optimization_settings.h" +#include "../helpers/make_time_grid.h" + +#include +#include +#include +#include +#include +#include + +template +void set_control_dampings(const SecirvvsOptimization& settings, mio::osecirvvs::Model& model, + const std::vector& control_parameter_values) +{ + std::vector time_steps = make_time_grid(settings.t0(), settings.tmax(), settings.num_intervals()); + + mio::UncertainContactMatrix& contacts = model.parameters.template get>(); + std::vector>& contact_dampings = contacts.get_dampings(); + + for (size_t control_interval = 0; control_interval < settings.num_control_intervals(); control_interval++) { + mio::SimulationTime time(time_steps[control_interval * settings.pc_resolution()]); + + for (size_t i = 0; i < settings.num_control_parameters(); i++) { + FP control_value = control_parameter_values[control_interval * settings.num_control_parameters() + i]; + + mio::DampingSampling damping = + settings.control_parameters()[i].make_damping_sampling(control_value, time); + + contact_dampings.push_back(damping); + } + } + contacts.make_matrix(); +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/create_constraints_ESID.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/create_constraints_ESID.h new file mode 100644 index 0000000000..a96d1d5a01 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/create_constraints_ESID.h @@ -0,0 +1,73 @@ +#pragma once + +#include "boost/filesystem.hpp" +#include +#include + +#include "constraints/constraints.h" + +#include +#include +#include +#include + +inline void load_constraints_from_file(const boost::filesystem::path& constraints_file, + std::vector& path_constraints, + std::vector& terminal_constraints) +{ + try { + boost::property_tree::ptree root; + boost::property_tree::read_json(constraints_file.string(), root); + + // Helper lambda for reading constraint arrays + auto load_constraints = [](const boost::property_tree::ptree& array) { + std::vector constraints; + for (const auto& item : array) { + const auto& obj = item.second; + + std::string name = obj.get("name"); + double min_val = obj.get("min"); + double max_val = obj.get("max"); + + // Handle optional or null node_id + std::optional node_id_opt = std::nullopt; + auto node_id_val = obj.get_optional("node_id"); + if (node_id_val.has_value()) { + node_id_opt = std::optional(*node_id_val); + } + + constraints.emplace_back(name, std::make_pair(min_val, max_val), node_id_opt); + } + return constraints; + }; + + // Load path_constraints and terminal_constraints (if present) + if (auto pc = root.get_child_optional("path_constraints")) { + path_constraints = load_constraints(*pc); + } + + if (auto tc = root.get_child_optional("terminal_constraints")) { + terminal_constraints = load_constraints(*tc); + } + + // Optional: Print debug info + std::cout << "Loaded " << path_constraints.size() << " path constraints.\n"; + for (const auto& c : path_constraints) { + std::cout << " - " << c.name() << " [" << c.min() << ", " << c.max() << "]"; + if (c.node_id()) + std::cout << " node_id=" << *c.node_id(); + std::cout << "\n"; + } + + std::cout << "Loaded " << terminal_constraints.size() << " terminal constraints.\n"; + for (const auto& c : terminal_constraints) { + std::cout << " - " << c.name() << " [" << c.min() << ", " << c.max() << "]"; + if (c.node_id()) + std::cout << " node_id=" << *c.node_id(); + std::cout << "\n"; + } + } + catch (const std::exception& e) { + std::cerr << "Error reading constraints from '" << constraints_file.string() << "': " << e.what() << std::endl; + } +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/create_control_parameters.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/create_control_parameters.h new file mode 100644 index 0000000000..7b8d87ab51 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/create_control_parameters.h @@ -0,0 +1,111 @@ +#pragma once + +#include "boost/filesystem.hpp" + +#include "memilio/epidemiology/damping.h" + +#include "control_parameters/control_parameters.h" + +#include +#include +#include + +enum class Intervention +{ + Home, + SchoolClosure, + HomeOffice, + GatheringBanFacilitiesClosure, + PhysicalDistanceAndMasks, + SeniorAwareness, + Count +}; + +enum class InterventionLevel +{ + Main, + PhysicalDistanceAndMasks, + SeniorAwareness, + Holidays, + Count +}; + +enum class ContactLocation +{ + Home, + School, + Work, + Other, + Count +}; + +std::vector create_control_parameters(const boost::filesystem::path& interventions_file, + size_t num_age_groups) +{ + std::vector control_parameters; + + try { + boost::property_tree::ptree pt; + boost::property_tree::read_json(interventions_file.string(), pt); + + std::cout << "Available interventions from '" << interventions_file.string() << "':\n"; + std::set intervention_names; + for (const auto& item : pt) { + try { + std::string name = item.second.get("name"); + intervention_names.insert(name); + std::cout << " - " << name << "\n"; + } + catch (...) { + // ignore malformed entries + } + } + + // Helper lambda to reduce duplication + auto add_control_parameter = [&](const std::string& name, const std::pair& allowed_range, + double effectiveness_value, double cost_value, InterventionLevel level_enum, + Intervention type_enum, const std::vector& locations_vec) { + if (intervention_names.find(name) != intervention_names.end()) { + mio::DampingLevel level(static_cast(level_enum)); + mio::DampingType type(static_cast(type_enum)); + + std::vector locations; + for (auto loc : locations_vec) { + locations.push_back(static_cast(loc)); + } + + Eigen::VectorX group_weights = Eigen::VectorX::Constant(num_age_groups, 1.0); + + control_parameters.emplace_back(name, allowed_range, effectiveness_value, cost_value, level, type, + locations, group_weights); + + std::cout << "Added control parameter for " << name << " intervention.\n"; + } + }; + + // clang-format off + // Add control parameters here + // name, range, effectiveness, cost, intervention level, intervention type, locations + add_control_parameter("School closure", {0.0, 1.0}, 1.0, 1.0, + InterventionLevel::Main, Intervention::SchoolClosure, {ContactLocation::School}); + + add_control_parameter("Face masks & social distancing School", {0.0, 1.0}, 0.25, 1.0, + InterventionLevel::PhysicalDistanceAndMasks, Intervention::PhysicalDistanceAndMasks, {ContactLocation::School}); + + add_control_parameter("Face masks & social distancing Work", {0.0, 1.0}, 0.25, 1.0, + InterventionLevel::PhysicalDistanceAndMasks, Intervention::PhysicalDistanceAndMasks, {ContactLocation::Work}); + + add_control_parameter("Face masks & social distancing Other", {0.0, 1.0}, 0.35, 1.0, + InterventionLevel::PhysicalDistanceAndMasks, Intervention::PhysicalDistanceAndMasks, {ContactLocation::Other}); + + add_control_parameter("Remote work", {0.0, 1.0}, 0.25, 1.0, + InterventionLevel::Main, Intervention::HomeOffice, {ContactLocation::Work}); + // clang-format on + } + catch (const std::exception& error) { + std::cerr << "Failed to read interventions file '" << interventions_file.string() << "': " << error.what() + << "\n"; + } + + return control_parameters; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/ad_type.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/ad_type.h new file mode 100644 index 0000000000..961ae6db99 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/ad_type.h @@ -0,0 +1,7 @@ +#pragma once + +enum class ADType +{ + Forward, + Reverse +}; diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/integrator_selector.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/integrator_selector.h new file mode 100644 index 0000000000..98e71f8f7f --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/integrator_selector.h @@ -0,0 +1,74 @@ +#pragma once + +#include "config.h" +#include "memilio/math/integrator.h" +#include "memilio/math/euler.h" +#include "memilio/math/adapt_rk.h" +#include "boost/numeric/odeint.hpp" + +enum class IntegratorType +{ + ExplicitEuler, + ExplicitCashKarp54, + ExplicitFehlberg78, + ControlledCashKarp54, + ControlledFehlberg78, + RK_Integrator +}; + +template +std::unique_ptr> make_integrator(IntegratorType type, FP dt) +{ + + FP dt_min = std::numeric_limits::min(); + FP dt_max = dt; + FP abs_tol = 1e-12; + FP rel_tol = 1e-8; + + using namespace boost::numeric::odeint; + + switch (type) { + case IntegratorType::ExplicitEuler: { + auto eulerIntegrator = std::make_unique>(); + return eulerIntegrator; + } + case IntegratorType::ExplicitCashKarp54: { + auto cashKarpStepper = std::make_unique>(); + cashKarpStepper->get_dt_min() = dt_min; + cashKarpStepper->get_dt_max() = dt_max; + return cashKarpStepper; + } + case IntegratorType::ExplicitFehlberg78: { + auto fehlbergStepper = std::make_unique>(); + fehlbergStepper->get_dt_min() = dt_min; + fehlbergStepper->get_dt_max() = dt_max; + return fehlbergStepper; + } + case IntegratorType::ControlledCashKarp54: { + auto cashKarpControlled = std::make_unique>(); + cashKarpControlled->set_dt_min(dt_min); + cashKarpControlled->set_dt_max(dt_max); + cashKarpControlled->set_abs_tolerance(abs_tol); + cashKarpControlled->set_rel_tolerance(rel_tol); + return cashKarpControlled; + } + case IntegratorType::ControlledFehlberg78: { + auto fehlbergControlled = std::make_unique>(); + fehlbergControlled->set_dt_min(dt_min); + fehlbergControlled->set_dt_max(dt_max); + fehlbergControlled->set_abs_tolerance(abs_tol); + fehlbergControlled->set_rel_tolerance(rel_tol); + return fehlbergControlled; + } + case IntegratorType::RK_Integrator: { + auto RK_Integrator = std::make_unique>(); + RK_Integrator->set_dt_min(dt_min); + RK_Integrator->set_dt_max(dt_max); + RK_Integrator->set_abs_tolerance(abs_tol); + RK_Integrator->set_rel_tolerance(rel_tol); + return RK_Integrator; + } + default: + throw std::invalid_argument("Unknown IntegratorType."); + } +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/make_time_grid.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/make_time_grid.h new file mode 100644 index 0000000000..322dbfc3d0 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/helpers/make_time_grid.h @@ -0,0 +1,15 @@ +#pragma once + +#include +#include + +template +std::vector make_time_grid(FP t0, FP tmax, size_t num_intervals) +{ + std::vector grid(num_intervals + 1); + FP grid_spacing = (tmax - t0) / num_intervals; + for (size_t i = 0; i <= num_intervals; i++) { + grid[i] = t0 + i * grid_spacing; + } + return grid; +} \ No newline at end of file diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/constraint_function.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/constraint_function.h new file mode 100644 index 0000000000..b8381a0630 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/constraint_function.h @@ -0,0 +1,71 @@ +#pragma once + +#include "../optimization_settings/optimization_settings.h" +#include "models/ode_secirvvs/model.h" + +#include "../control_parameters/damping_controls.h" +#include "../constraints/update_constraints.h" + +#include "../helpers/integrator_selector.h" +#include "../helpers/make_time_grid.h" + +#include "../constraints/infection_state_utils.h" + +#include + +template +void constraint_function( + mio::Graph>, mio::MobilityEdge> graph_model, + const SecirvvsOptimization& settings, const FP* ptr_parameters, size_t n, FP* ptr_constraints, size_t m) +{ + // ----------------------------------------------------------------- // + // Evaluate the constraints on the model. // + // Step 1. Define dampings based on 'FP* ptr_parameters'. // + // Step 2. Gather information in 'path_values' in 'terminal_values'. // + // Step 3. Fill information into 'FP* ptr_constraints'. // + // ----------------------------------------------------------------- // + assert(n == settings.num_control_parameters() * settings.num_control_intervals()); + assert(m == settings.num_path_constraints() + settings.num_terminal_constraints()); + + std::vector path_constraint_values(settings.num_path_constraints(), 0.0); + std::vector terminal_constraint_values(settings.num_terminal_constraints(), 0.0); + + std::vector parameters(settings.num_control_parameters() * settings.num_control_intervals()); + for (size_t control_interval = 0; control_interval < settings.num_control_intervals(); control_interval++) { + for (size_t control_index = 0; control_index < settings.num_control_parameters(); control_index++) { + size_t idx = control_interval * settings.num_control_parameters() + control_index; + parameters[idx] = settings.activation_function()(ptr_parameters[idx]); + } + } + + for (auto& node : graph_model.nodes()) { + set_control_dampings(settings, node.property.get_simulation().get_model(), parameters); + } + + for (auto& node : graph_model.nodes()) { + node.property.get_simulation().set_integrator_core( + std::move(make_integrator(settings.integrator_type(), settings.dt()))); + } + + std::vector time_steps = make_time_grid(settings.t0(), settings.tmax(), settings.num_intervals()); + auto graph_simulation = mio::make_mobility_sim(settings.t0(), settings.dt(), graph_model); + + update_path_constraint(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(settings, graph_simulation.get_graph(), path_constraint_values); + } + update_terminal_constraint(settings, graph_simulation.get_graph(), terminal_constraint_values); + + size_t idx = 0; + for (size_t path_constraint_index = 0; path_constraint_index < settings.num_path_constraints(); + path_constraint_index++) { + ptr_constraints[idx] = path_constraint_values[path_constraint_index]; + idx++; + } + for (size_t terminal_constraint_index = 0; terminal_constraint_index < settings.num_terminal_constraints(); + terminal_constraint_index++) { + ptr_constraints[idx] = terminal_constraint_values[terminal_constraint_index]; + idx++; + } +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/objective_function.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/objective_function.h new file mode 100644 index 0000000000..ebd1ae86a4 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/objective_function.h @@ -0,0 +1,68 @@ +#pragma once + +#include "../optimization_settings/optimization_settings.h" +#include "models/ode_secirvvs/model.h" + +#include "../control_parameters/damping_controls.h" +#include "../constraints/update_constraints.h" + +#include "../helpers/integrator_selector.h" +#include "../helpers/make_time_grid.h" + +#include "../constraints/infection_state_utils.h" + +#include + +template +FP objective_function( + mio::Graph>, mio::MobilityEdge> graph_model, + const SecirvvsOptimization& settings, const FP* ptr_parameters, size_t n) +{ + // ------------------------------------------------------------------ // + // Evaluate the objective function of the model. // + // Step 1. Define dampings based on 'const FP* parameters'. // + // Step 2. Evaluate the objective function based on // + // the parameters and the infection states in the simulation. // + // ------------------------------------------------------------------ // + assert(n == settings.num_control_parameters() * settings.num_control_intervals()); + + FP objective = 0.0; + + std::vector parameters(settings.num_control_parameters() * settings.num_control_intervals()); + + for (size_t control_interval = 0; control_interval < settings.num_control_intervals(); control_interval++) { + for (size_t control_index = 0; control_index < settings.num_control_parameters(); control_index++) { + size_t idx = control_interval * settings.num_control_parameters() + control_index; + parameters[idx] = settings.activation_function()(ptr_parameters[idx]); + } + } + + for (auto& node : graph_model.nodes()) { + set_control_dampings(settings, node.property.get_simulation().get_model(), parameters); + } + + for (auto& node : graph_model.nodes()) { + node.property.get_simulation().set_integrator_core( + std::move(make_integrator(settings.integrator_type(), settings.dt()))); + } + + std::vector time_steps = make_time_grid(settings.t0(), settings.tmax(), settings.num_intervals()); + auto graph_simulation = mio::make_mobility_sim(settings.t0(), settings.dt(), graph_model); + + for (size_t interval = 0; interval < settings.num_intervals(); interval++) { + + graph_simulation.advance(time_steps[interval + 1]); + + size_t control_interval = interval / settings.pc_resolution(); + + // Add costs of all control parameters for current interval + FP interval_cost = 0.0; + for (size_t control_index = 0; control_index < settings.num_control_parameters(); control_index++) { + interval_cost += settings.control_parameters()[control_index].cost() * + parameters[control_interval * settings.num_control_parameters() + control_index]; + } + objective += interval_cost / settings.num_intervals(); + } + + return objective; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/save_solution.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/save_solution.h new file mode 100644 index 0000000000..0a79ecba24 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/save_solution.h @@ -0,0 +1,85 @@ + +#pragma once + +#include "../optimization_settings/optimization_settings.h" +#include "models/ode_secirvvs/model.h" + +#include "../control_parameters/damping_controls.h" +#include "../constraints/update_constraints.h" + +#include "../helpers/integrator_selector.h" +#include "../helpers/make_time_grid.h" + +#include "../constraints/infection_state_utils.h" + +#include +#include + +template +void save_solution( + mio::Graph>, mio::MobilityEdge> graph_model, + const SecirvvsOptimization& settings, size_t n, const FP* ptr_parameters, const FP* z_L, const FP* z_U, size_t m, + const FP* ptr_constraints, const FP* lambda, FP obj_value) +{ + std::vector parameters(settings.num_control_parameters() * settings.num_control_intervals()); + + for (size_t control_interval = 0; control_interval < settings.num_control_intervals(); control_interval++) { + for (size_t control_index = 0; control_index < settings.num_control_parameters(); control_index++) { + size_t idx = control_interval * settings.num_control_parameters() + control_index; + + parameters[idx] = settings.activation_function()(ptr_parameters[idx]); + } + } + + std::cout << "\nPath Constraints:\n"; + for (size_t i = 0; i < settings.num_path_constraints(); ++i) { + std::cout << settings.path_constraints()[i].name() << ": " << ptr_constraints[i] << std::endl; + } + std::cout << "\nTerminal Constraints:\n"; + for (size_t i = 0; i < settings.num_terminal_constraints(); ++i) { + std::cout << settings.terminal_constraints()[i].name() << ": " + << ptr_constraints[i + settings.num_path_constraints()] << std::endl; + } + + // Write parameters to CSV file + std::ofstream file("control_parameters.csv"); + if (!file.is_open()) { + std::cerr << "Failed to open control_parameters.csv for writing.\n"; + return; + } + + // Write header + file << "Time"; + for (const auto& cp : settings.control_parameters()) { + file << ", " << cp.name(); + } + file << "\n"; + + size_t num_intervals = settings.num_control_intervals(); + size_t num_controls = settings.num_control_parameters(); + double dt = (settings.tmax() - settings.t0()) / num_intervals; + + // Write data rows + for (size_t i = 0; i < num_intervals; ++i) { + double time = settings.t0() + i * dt; + file << time; + + size_t base_index = i * num_controls; + for (size_t j = 0; j < num_controls; ++j) { + file << ", " << parameters[base_index + j]; + } + file << "\n"; + } + + // Repeat last set of parameters at final time + double final_time = settings.t0() + num_intervals * dt; + file << final_time; + + size_t last_base_index = (num_intervals - 1) * num_controls; + for (size_t j = 0; j < num_controls; ++j) { + file << ", " << parameters[last_base_index + j]; + } + file << "\n"; + + file.close(); +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/secirvvs_ipopt.cpp b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/secirvvs_ipopt.cpp new file mode 100644 index 0000000000..7ce9836931 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/secirvvs_ipopt.cpp @@ -0,0 +1,403 @@ +#include "secirvvs_ipopt.h" + +#include "IpIpoptData.hpp" + +#include "memilio/utils/mioomp.h" + +#include "objective_function.h" +#include "constraint_function.h" +#include "save_solution.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +Secirvvs_NLP::Secirvvs_NLP(const SecirvvsOptimization& settings) + : m_settings(settings) +{ + size_t num_control_intervals = m_settings.num_control_intervals(); + size_t num_control_parameters = m_settings.num_control_parameters(); + size_t num_path_constraints = m_settings.num_path_constraints(); + size_t num_terminal_constraints = m_settings.num_terminal_constraints(); + + m_n = num_control_intervals * num_control_parameters; + m_m = num_path_constraints + num_terminal_constraints; + + m_nnz_jac_g = m_n * m_m; + bool use_exact_hessian = false; + m_nnz_h_lag = use_exact_hessian ? m_n * (m_n + 1) / 2 : 0; + + // Resize vectors + m_x_l.resize(m_n); + m_x_u.resize(m_n); + m_g_l.resize(m_m); + m_g_u.resize(m_m); + + // Fill m_x_l and m_x_u from control parameters + for (size_t control_interval = 0; control_interval < num_control_intervals; control_interval++) { + for (size_t control_index = 0; control_index < num_control_parameters; control_index++) { + size_t idx = control_interval * num_control_parameters + control_index; + m_x_l[idx] = m_settings.control_parameters()[control_index].min(); + m_x_u[idx] = m_settings.control_parameters()[control_index].max(); + } + } + // Fill m_g_l and m_g_u from path- and terminal constraints + size_t constraint_index = 0; + for (size_t path_constraint_index = 0; path_constraint_index < num_path_constraints; path_constraint_index++) { + m_g_l[constraint_index] = m_settings.path_constraints()[path_constraint_index].min(); + m_g_u[constraint_index] = m_settings.path_constraints()[path_constraint_index].max(); + constraint_index++; + } + for (size_t terminal_constraint_index = 0; terminal_constraint_index < num_terminal_constraints; + terminal_constraint_index++) { + m_g_l[constraint_index] = m_settings.terminal_constraints()[terminal_constraint_index].min(); + m_g_u[constraint_index] = m_settings.terminal_constraints()[terminal_constraint_index].max(); + constraint_index++; + } + + if (m_settings.ad_eval_f() == ADType::Reverse || m_settings.ad_eval_jac() == ADType::Reverse) { + if (!ad::ga1s::global_tape) { + ad::ga1s::global_tape = tape_t::create(); + } + m_tape = ad::ga1s::global_tape; + } +} + +Secirvvs_NLP::~Secirvvs_NLP() +{ + if (m_tape) { + tape_t::remove(m_tape); + m_tape = nullptr; + } +} + +bool Secirvvs_NLP::get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g, Ipopt::Index& nnz_h_lag, + Ipopt::TNLP::IndexStyleEnum& index_style) +{ + n = m_n; + m = m_m; + nnz_jac_g = m_nnz_jac_g; + nnz_h_lag = m_nnz_h_lag; + index_style = Ipopt::TNLP::C_STYLE; + return true; +} + +bool Secirvvs_NLP::get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u, Ipopt::Index m, + Ipopt::Number* g_l, Ipopt::Number* g_u) +{ + assert(n == m_n && m == m_m); + + for (Ipopt::Index i = 0; i < n; i++) { + x_l[i] = m_x_l[i]; + x_u[i] = m_x_u[i]; + } + + for (Ipopt::Index i = 0; i < m; i++) { + g_l[i] = m_g_l[i]; + g_u[i] = m_g_u[i]; + } + + return true; +} + +bool Secirvvs_NLP::get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x, bool init_z, + Ipopt::Number* /*z_L*/, Ipopt::Number* /*z_U*/, Ipopt::Index m, bool init_lambda, + Ipopt::Number* /*lambda*/ +) +{ + assert(init_x && !init_z && !init_lambda); + assert(n == m_n && m == m_m); + + mio::unused(n, m, init_x, init_z, init_lambda); + + size_t num_control_intervals = m_settings.num_control_intervals(); + size_t num_control_parameters = m_settings.num_control_parameters(); + + std::random_device rd; + std::mt19937 gen(rd()); + + for (size_t control_interval = 0; control_interval < num_control_intervals; control_interval++) { + for (size_t control_index = 0; control_index < num_control_parameters; control_index++) { + size_t idx = control_interval * num_control_parameters + control_index; + + double min_val = m_settings.control_parameters()[control_index].min(); + double max_val = m_settings.control_parameters()[control_index].max(); + + if (m_settings.random_start()) { + std::uniform_real_distribution<> dis(min_val, max_val); + x[idx] = dis(gen); + } + else { + x[idx] = 0.5 * (min_val + max_val); + } + } + } + + return true; +} + +bool Secirvvs_NLP::eval_f(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Number& obj_value) +{ + assert(n == m_n); + + auto graph_model = m_settings.optimization_model().get_graph_model(); + obj_value = objective_function(graph_model, m_settings, x, n); + + return true; +} + +bool Secirvvs_NLP::eval_g(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Index m, Ipopt::Number* g) +{ + assert(n == m_n && m == m_m); + + auto graph_model = m_settings.optimization_model().get_graph_model(); + constraint_function(graph_model, m_settings, x, n, g, m); + + return true; +} + +bool Secirvvs_NLP::eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f) +{ + assert(n == m_n); + + if (m_settings.ad_eval_f() == ADType::Forward) { + // Use forward mode for Jacobian evaluation + return eval_grad_f_forward(n, x, new_x, grad_f); + } + else { + // Use reverse mode for Jacobian evaluation + return eval_grad_f_reverse(n, x, new_x, grad_f); + } +} + +bool Secirvvs_NLP::eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Index nele_jac, + Ipopt::Index* iRow, Ipopt::Index* jCol, Ipopt::Number* values) +{ + assert(n == m_n && m == m_m && nele_jac == m_nnz_jac_g); + + if (m_settings.ad_eval_jac() == ADType::Forward) { + // Use forward mode for Jacobian evaluation + return eval_jac_g_forward(n, x, new_x, m, nele_jac, iRow, jCol, values); + } + else { + // Use reverse mode for Jacobian evaluation + return eval_jac_g_reverse(n, x, new_x, m, nele_jac, iRow, jCol, values); + } +} + +bool Secirvvs_NLP::eval_h(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Number obj_factor, + Ipopt::Index m, const Ipopt::Number* lambda, bool /*new_lambda*/, Ipopt::Index nele_hess, + Ipopt::Index* iRow, Ipopt::Index* jCol, Ipopt::Number* values) +{ + return false; +} + +bool Secirvvs_NLP::eval_grad_f_forward(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Number* grad_f) +{ + assert(n == m_n); + using ad_t = ad::gt1s::type; + + // Define graph model outside of parallel region + auto graph_model = m_settings.optimization_model().get_graph_model(); + + PRAGMA_OMP(parallel) + { + std::vector x_ad(n); + for (Ipopt::Index i = 0; i < n; ++i) { + x_ad[i] = x[i]; + ad::derivative(x_ad[i]) = 0.0; + } + + PRAGMA_OMP(for) + for (Ipopt::Index column = 0; column < n; ++column) { + ad::derivative(x_ad[column]) = 1.0; + ad_t obj_ad = objective_function(graph_model, m_settings, x_ad.data(), n); + grad_f[column] = ad::derivative(obj_ad); + ad::derivative(x_ad[column]) = 0.0; + } + } + + return true; +} + +bool Secirvvs_NLP::eval_grad_f_reverse(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Number* grad_f) +{ + assert(n == m_n); + using ad_t = ad::ga1s::type; + + auto graph_model = m_settings.optimization_model().get_graph_model(); + + m_tape->reset(); + + std::vector x_ad(n); + + for (Ipopt::Index i = 0; i < n; ++i) { + ad::value(x_ad[i]) = x[i]; + ad::derivative(x_ad[i]) = 0.0; + m_tape->register_variable(x_ad[i]); + } + + ad_t obj_ad = objective_function(graph_model, m_settings, x_ad.data(), n); + m_tape->register_output_variable(obj_ad); + ad::derivative(obj_ad) = 1.0; + + m_tape->interpret_adjoint(); + + for (Ipopt::Index i = 0; i < n; ++i) { + grad_f[i] = ad::derivative(x_ad[i]); + } + + return true; +} + +bool Secirvvs_NLP::eval_jac_g_forward(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Index m, + Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index* jCol, + Ipopt::Number* values) +{ + assert(n == m_n && m == m_m && nele_jac == m_nnz_jac_g); + + if (values == nullptr) { + Ipopt::Index idx = 0; + for (Ipopt::Index column = 0; column < n; ++column) { + for (Ipopt::Index row = 0; row < m; ++row) { + iRow[idx] = row; + jCol[idx] = column; + ++idx; + } + } + return true; + } + + using ad_t = ad::gt1s::type; + + // Define graph model outside of parallel region + auto graph_model = m_settings.optimization_model().get_graph_model(); + + PRAGMA_OMP(parallel) + { + std::vector x_ad(n); + std::vector g_ad(m); + + for (Ipopt::Index i = 0; i < n; ++i) { + x_ad[i] = x[i]; + ad::derivative(x_ad[i]) = 0.0; + } + + PRAGMA_OMP(for) + for (Ipopt::Index column = 0; column < n; ++column) { + ad::derivative(x_ad[column]) = 1.0; + constraint_function(graph_model, m_settings, x_ad.data(), n, g_ad.data(), m); + + for (Ipopt::Index row = 0; row < m; ++row) { + values[column * m + row] = ad::derivative(g_ad[row]); + } + + ad::derivative(x_ad[column]) = 0.0; + } + } + + // Critical check: Abort if gradient is too large + const double gradient_threshold = 1e15; + for (Ipopt::Index column = 0; column < n; ++column) { + for (Ipopt::Index row = 0; row < m; ++row) { + const double grad_val = values[column * m + row]; + if (std::abs(grad_val) > gradient_threshold) { + mio::log_critical("Gradient explosion detected: |∂g_{}/∂x_{}| = {:e} at x[{}] = {:e}", row, column, + grad_val, column, x[column]); + return false; + } + } + } + + return true; +} + +bool Secirvvs_NLP::eval_jac_g_reverse(Ipopt::Index n, const Ipopt::Number* x, bool /*new_x*/, Ipopt::Index m, + Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index* jCol, + Ipopt::Number* values) +{ + assert(n == m_n && m == m_m && nele_jac == m_nnz_jac_g); + + if (values == nullptr) { + Ipopt::Index idx = 0; + for (Ipopt::Index row = 0; row < m; ++row) { + for (Ipopt::Index column = 0; column < n; ++column) { + iRow[idx] = row; + jCol[idx] = column; + ++idx; + } + } + return true; + } + + using ad_t = ad::ga1s::type; + + auto graph_model = m_settings.optimization_model().get_graph_model(); + + m_tape->reset(); + + std::vector x_ad(n); + for (Ipopt::Index column = 0; column < n; ++column) { + ad::value(x_ad[column]) = x[column]; + ad::derivative(x_ad[column]) = 0.0; + m_tape->register_variable(x_ad[column]); + } + + std::vector g_ad(m); + constraint_function(graph_model, m_settings, x_ad.data(), n, g_ad.data(), m); + + for (Ipopt::Index row = 0; row < m; ++row) { + m_tape->register_output_variable(g_ad[row]); + } + + for (Ipopt::Index row = 0; row < m; ++row) { + m_tape->zero_adjoints(); + ad::derivative(g_ad[row]) = 1.0; + m_tape->interpret_adjoint(); // Takes up most of the time + for (Ipopt::Index column = 0; column < n; ++column) { + values[row * n + column] = ad::derivative(x_ad[column]); + } + } + + // Critical check: Abort if gradient is too large + const double gradient_threshold = 1e15; + for (Ipopt::Index row = 0; row < m; ++row) { + for (Ipopt::Index column = 0; column < n; ++column) { + const double grad_val = values[row * n + column]; + if (std::abs(grad_val) > gradient_threshold) { + mio::log_critical("Gradient explosion detected: |∂g_{}/∂x_{}| = {:e} at x[{}] = {:e}", row, column, + grad_val, column, x[column]); + return false; + } + } + } + + return true; +} + +void Secirvvs_NLP::finalize_solution(Ipopt::SolverReturn status, Ipopt::Index n, const Ipopt::Number* x, + const Ipopt::Number* z_L, const Ipopt::Number* z_U, Ipopt::Index m, + const Ipopt::Number* g, const Ipopt::Number* lambda, Ipopt::Number obj_value, + const Ipopt::IpoptData* ip_data, Ipopt::IpoptCalculatedQuantities* ip_cq) +{ + std::cout << "Final Objective Value: " << obj_value << "\n"; + + if (ip_data) { + std::cout << "Number of iterations: " << ip_data->iter_count() << "\n"; + } + + auto graph_model = m_settings.optimization_model().get_graph_model(); + save_solution(graph_model, m_settings, n, x, z_L, z_U, m, g, lambda, obj_value); + + return; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/secirvvs_ipopt.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/secirvvs_ipopt.h new file mode 100644 index 0000000000..9d2b1458e6 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/ipopt_solver/secirvvs_ipopt.h @@ -0,0 +1,71 @@ +#pragma once + +#include "IpTNLP.hpp" +#include "memilio/ad/ad.h" + +#include "../optimization_settings/optimization_settings.h" +#include "models/ode_secirvvs/model.h" + +#include + +class Secirvvs_NLP : public Ipopt::TNLP +{ +public: + Secirvvs_NLP(const SecirvvsOptimization& settings); + + Secirvvs_NLP(const Secirvvs_NLP&) = delete; + Secirvvs_NLP(Secirvvs_NLP&&) = delete; + Secirvvs_NLP& operator=(const Secirvvs_NLP&) = delete; + Secirvvs_NLP& operator=(Secirvvs_NLP&&) = delete; + + ~Secirvvs_NLP() override; + + bool get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g, Ipopt::Index& nnz_h_lag, + Ipopt::TNLP::IndexStyleEnum& index_style) override; + + bool get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u, Ipopt::Index m, Ipopt::Number* g_l, + Ipopt::Number* g_u) override; + + bool get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x, bool init_z, Ipopt::Number* z_L, + Ipopt::Number* z_U, Ipopt::Index m, bool init_lambda, Ipopt::Number* lambda) override; + + bool eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value) override; + bool eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f) override; + + bool eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g) override; + bool eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Index nele_jac, + Ipopt::Index* iRow, Ipopt::Index* jCol, Ipopt::Number* values) override; + + bool eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number obj_factor, Ipopt::Index m, + const Ipopt::Number* lambda, bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow, + Ipopt::Index* jCol, Ipopt::Number* values) override; + + bool eval_grad_f_forward(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f); + bool eval_grad_f_reverse(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f); + + bool eval_jac_g_forward(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Index nele_jac, + Ipopt::Index* iRow, Ipopt::Index* jCol, Ipopt::Number* values); + bool eval_jac_g_reverse(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Index nele_jac, + Ipopt::Index* iRow, Ipopt::Index* jCol, Ipopt::Number* values); + + void finalize_solution(Ipopt::SolverReturn status, Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, + const Ipopt::Number* z_U, Ipopt::Index m, const Ipopt::Number* g, + const Ipopt::Number* lambda, Ipopt::Number obj_value, const Ipopt::IpoptData* ip_data, + Ipopt::IpoptCalculatedQuantities* ip_cq) override; + +private: + const SecirvvsOptimization& m_settings; + + Ipopt::Index m_n; + Ipopt::Index m_m; + Ipopt::Index m_nnz_jac_g; + Ipopt::Index m_nnz_h_lag; + + std::vector m_x_l; + std::vector m_x_u; + std::vector m_g_l; + std::vector m_g_u; + + using tape_t = ad::ga1s::tape_t; + tape_t* m_tape = nullptr; +}; diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_ESID.cpp b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_ESID.cpp new file mode 100644 index 0000000000..91bd7a0adb --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_ESID.cpp @@ -0,0 +1,136 @@ +#include "boost/filesystem.hpp" +#include +#include + +#include "memilio/utils/logging.h" +#include "memilio/math/eigen.h" +#include "memilio/epidemiology/damping.h" + +#include "create_control_parameters.h" +#include "create_constraints_ESID.h" +#include "check_feasibility.h" + +#include "optimization_model/optimization_model.h" +#include "optimization_settings/optimization_settings.h" + +#include "helpers/integrator_selector.h" +#include "helpers/ad_type.h" +#include "control_parameters/control_activation.h" +#include "constraints/constraints.h" + +#include "ipopt_solver/secirvvs_ipopt.h" +#include "IpTNLP.hpp" +#include "IpIpoptApplication.hpp" + +#include +#include +#include +#include + +int main(int argc, char* argv[]) +{ + if (argc != 3) { + std::cerr << "Usage: " << argv[0] << " [data directory] [optimal_control_settings directory]\n"; + return 1; + } + + mio::set_log_level(mio::LogLevel::critical); + + boost::filesystem::path data_directory = argv[1]; + boost::filesystem::path optimal_control_settings_directory = argv[2]; + + boost::filesystem::path interventions_file = optimal_control_settings_directory / "intervention_list.json"; + boost::filesystem::path constraints_file = optimal_control_settings_directory / "constraints.json"; + + size_t num_age_groups = 6; + std::vector control_parameters = + create_control_parameters(interventions_file, num_age_groups); + + std::vector path_constraints; + std::vector terminal_constraints; + load_constraints_from_file(constraints_file, path_constraints, terminal_constraints); + + size_t simulation_weeks = 8; + size_t days_per_week = 7; + size_t simulation_days = simulation_weeks * days_per_week; + + OptimizationModel optimization_model(data_directory, simulation_days, num_age_groups); + + auto graph_model = optimization_model.get_graph_model(); + const auto& nodes = graph_model.nodes(); + + std::cout << "Number of graph nodes: " << nodes.size() << std::endl; + + // --------------------------------------- // + // --- Configure optimization settings --- // + // --------------------------------------- // + size_t num_control_intervals = simulation_weeks; // Number of control intervals (e.g., one for each week) + size_t pc_resolution = days_per_week; // Path constraint resolution (e.g., 7 days per control interval) + bool random_start = false; // Set to true to randomize initial control values + + IntegratorType integrator_type = IntegratorType::ControlledCashKarp54; + size_t integrator_resolution = 10; + + // ADType::Reverse can run out of tape size: Segmentation fault + ADType ad_eval_f = ADType::Forward; + ADType ad_eval_jac = ADType::Forward; + + // Mapping f:[0,1]->[0,1], f(0)=0, f(1/2)=1/2, f(1)=1. + // ControlActivation::Sigmoid can only be used if controls are in [0,1]. + ControlActivation activation = ControlActivation::Linear; + + SecirvvsOptimization settings(optimization_model, num_control_intervals, pc_resolution, random_start, + integrator_type, integrator_resolution, ad_eval_f, ad_eval_jac, control_parameters, + path_constraints, terminal_constraints, activation); + + check_constraint_feasibility(settings); + + // ----------------------------- // + // --- Create NLP and solver --- // + // ----------------------------- // + Ipopt::SmartPtr mynlp = new Secirvvs_NLP(settings); + Ipopt::SmartPtr app = IpoptApplicationFactory(); + + int verbose = 5; + int print_frequency = 1; + bool print_timings = true; + int max_iter = 500; + double tol = 1e-5; + bool use_exact_hessian = false; + assert(!use_exact_hessian); + std::string hessianApproximation = use_exact_hessian ? "exact" : "limited-memory"; + std::string linearSolver = "mumps"; + std::string muStrategy = "adaptive"; + std::string memoryUpdateType = "bfgs"; + + app->Options()->SetIntegerValue("print_level", verbose); + app->Options()->SetIntegerValue("print_frequency_iter", print_frequency); + app->Options()->SetIntegerValue("max_iter", max_iter); + app->Options()->SetNumericValue("tol", tol); + app->Options()->SetStringValue("linear_solver", linearSolver); + app->Options()->SetStringValue("hessian_approximation", hessianApproximation); + app->Options()->SetStringValue("mu_strategy", muStrategy); + app->Options()->SetStringValue("limited_memory_update_type", memoryUpdateType); + if (print_timings) { + app->Options()->SetStringValue("timing_statistics", "yes"); + app->Options()->SetStringValue("print_timing_statistics", "yes"); + } + + // Initialize and solve + Ipopt::ApplicationReturnStatus status = app->Initialize(); + if (status != Ipopt::Solve_Succeeded) { + std::cout << "\n*** Error during initialization!\n"; + return (int)status; + } + + status = app->OptimizeTNLP(mynlp); + + if (status == Ipopt::Solve_Succeeded) { + std::cout << "\n*** The problem was solved successfully!\n"; + } + else { + std::cout << "\n*** The problem failed to solve!\n"; + } + + return 0; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_settings/constraints.json b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_settings/constraints.json new file mode 100644 index 0000000000..f94f4e37ed --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_settings/constraints.json @@ -0,0 +1,30 @@ +{ + "path_constraints": [ + { + "name": "Infected", + "min": 0.0, + "max": 1000000.3, + "node_id": null + }, + { + "name": "Severe", + "min": 0.0, + "max": 125000.0, + "node_id": null + }, + { + "name": "Critical", + "min": 0.0, + "max": 200.0, + "node_id": null + } + ], + "terminal_constraints": [ + { + "name": "Dead", + "min": 0.0, + "max": 400000.0, + "node_id": null + } + ] +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_settings/intervention_list.json b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_settings/intervention_list.json new file mode 100644 index 0000000000..16488fc289 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimal_control_settings/intervention_list.json @@ -0,0 +1,32 @@ +[ + { + "id": "30513de3-e03a-469a-a3e3-635d45d1b2a4", + "name": "School closure", + "description": "School closure intervention", + "tags": null + }, + { + "id": "fb4fe730-5a7c-46a2-b731-494cfbd916b4", + "name": "Face masks & social distancing School", + "description": "Face mask usage and social distancing measures applied at 1-25% in schools", + "tags": null + }, + { + "id": "2dfb8d4e-cdbe-45ac-ae9c-d070fe8ece2e", + "name": "Face masks & social distancing Work", + "description": "Face mask usage and social distancing measures applied at 1-25% in workplaces", + "tags": null + }, + { + "id": "b5205110-96d3-4e0c-993c-2ca2dfe1412e", + "name": "Face masks & social distancing Other", + "description": "Face mask usage and social distancing measures applied at 1-25% in other settings", + "tags": null + }, + { + "id": "96cd1086-0c57-4e87-9977-bc43b2ba56f2", + "name": "Remote work", + "description": "Implementation of remote work policies", + "tags": null + } +] diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_model/optimization_model.cpp b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_model/optimization_model.cpp new file mode 100644 index 0000000000..b795ead29a --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_model/optimization_model.cpp @@ -0,0 +1,25 @@ +#include "optimization_model.h" + +OptimizationModel::OptimizationModel(const boost::filesystem::path& data_directory, size_t simulation_days, + size_t num_age_groups) + : m_data_directory(data_directory) + , m_simulation_days(simulation_days) + , m_num_age_groups(num_age_groups) + , m_cache(std::make_shared()) +{ + std::cout << "OptimizationModel initialized with data directory: " << m_data_directory << std::endl; +} + +boost::filesystem::path OptimizationModel::data_directory() const +{ + return this->m_data_directory; +} + +size_t OptimizationModel::simulation_days() const +{ + return this->m_simulation_days; +} +size_t OptimizationModel::num_age_groups() const +{ + return this->m_num_age_groups; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_model/optimization_model.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_model/optimization_model.h new file mode 100644 index 0000000000..fff99f0a68 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_model/optimization_model.h @@ -0,0 +1,292 @@ +#pragma once + +#include "boost/filesystem.hpp" + +#include "memilio/mobility/graph.h" +#include "memilio/mobility/metapopulation_mobility_instant.h" +#include "models/ode_secirvvs/model.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class OptimizationModel +{ +public: + OptimizationModel(const boost::filesystem::path& data_directory, size_t simulation_days, size_t num_age_groups); + + boost::filesystem::path data_directory() const; + size_t simulation_days() const; + size_t num_age_groups() const; + + // Retrieves (or generates) a graph model of the simulation. + // Uses template-based caching since reading in models can be time-consuming. + template + mio::Graph>, mio::MobilityEdge> get_graph_model() const; + +private: + boost::filesystem::path m_data_directory; + size_t m_simulation_days; + size_t m_num_age_groups; + + // Creates an artificial (synthetic) graph model with preset parameters and population data. + template + mio::Graph>, mio::MobilityEdge> + create_artificial_graph_model() const; + + struct Cache { + std::unordered_map data; + std::mutex mutex; + }; + std::shared_ptr m_cache; +}; + +// Retrieves a graph model from the cache or builds a new one if necessary. +template +auto OptimizationModel::get_graph_model() const + -> mio::Graph>, mio::MobilityEdge> +{ + const auto key = std::type_index(typeid(FP)); + { + // Check if cached + std::scoped_lock lock(m_cache->mutex); + auto it = m_cache->data.find(key); + if (it != m_cache->data.end()) { + return std::any_cast< + mio::Graph>, mio::MobilityEdge>>(it->second); + } + } + auto graph = create_artificial_graph_model(); + { + // Cache the result + std::scoped_lock lock(m_cache->mutex); + m_cache->data.emplace(key, graph); + } + return graph; +} + +// ------------------------------------------------------------ // +// Creates a synthetic graph model for testing the optimization // +template +mio::Graph>, mio::MobilityEdge> +OptimizationModel::create_artificial_graph_model() const +{ + // --------------------------- // + // Define the model parameters // + constexpr size_t num_age_groups = 6; + mio::osecirvvs::Model model(num_age_groups); + auto& params = model.parameters; + + params.template get>() = 100; + params.template get>() = 0.0143; + params.template get>() = 0.2; + params.template get>() = 7; + + params.template get>().resize( + mio::SimulationDay(m_simulation_days + 1)); + params.template get>().resize(mio::SimulationDay(m_simulation_days + 1)); + + size_t daily_vaccinations = 10; + for (size_t age_group = 0; age_group < num_age_groups; age_group++) { + for (size_t current_day = 0; current_day <= m_simulation_days; current_day++) { + mio::Index index{mio::AgeGroup(age_group), + mio::SimulationDay(current_day)}; + FP num_vaccinations = static_cast(current_day * daily_vaccinations); + params.template get>()[index] = num_vaccinations; + params.template get>()[index] = num_vaccinations; + } + } + + auto set_all_groups = [&](auto Tag, FP value) { + auto& age_group_params = params.template get(); + for (size_t age_group = 0; age_group < num_age_groups; age_group++) { + age_group_params[mio::AgeGroup(age_group)] = value; + } + }; + + // — times (all groups same value) + set_all_groups(mio::osecirvvs::TimeExposed{}, 3.33); + set_all_groups(mio::osecirvvs::TimeInfectedNoSymptoms{}, 1.87); + set_all_groups(mio::osecirvvs::TimeInfectedSymptoms{}, 7); + set_all_groups(mio::osecirvvs::TimeInfectedSevere{}, 6); + set_all_groups(mio::osecirvvs::TimeInfectedCritical{}, 7); + // — probabilities (all groups same value) + set_all_groups(mio::osecirvvs::TransmissionProbabilityOnContact{}, 0.15); + set_all_groups(mio::osecirvvs::RelativeTransmissionNoSymptoms{}, 0.5); + set_all_groups(mio::osecirvvs::RiskOfInfectionFromSymptomatic{}, 0.0); + set_all_groups(mio::osecirvvs::MaxRiskOfInfectionFromSymptomatic{}, 0.4); + set_all_groups(mio::osecirvvs::RecoveredPerInfectedNoSymptoms{}, 0.2); + set_all_groups(mio::osecirvvs::SeverePerInfectedSymptoms{}, 0.1); + set_all_groups(mio::osecirvvs::CriticalPerSevere{}, 0.1); + set_all_groups(mio::osecirvvs::DeathsPerCritical{}, 0.1); + // — immunity reductions (all groups same value) + set_all_groups(mio::osecirvvs::ReducExposedPartialImmunity{}, 0.8); + set_all_groups(mio::osecirvvs::ReducExposedImprovedImmunity{}, 0.331); + set_all_groups(mio::osecirvvs::ReducInfectedSymptomsPartialImmunity{}, 0.65); + set_all_groups(mio::osecirvvs::ReducInfectedSymptomsImprovedImmunity{}, 0.243); + set_all_groups(mio::osecirvvs::ReducInfectedSevereCriticalDeadPartialImmunity{}, 0.1); + set_all_groups(mio::osecirvvs::ReducInfectedSevereCriticalDeadImprovedImmunity{}, 0.091); + set_all_groups(mio::osecirvvs::ReducTimeInfectedMild{}, 0.9); + + // --------------------------- // + // Define the contact matrices // + int num_contact_locations = 4; + Eigen::Matrix baseline_home; + Eigen::Matrix baseline_school_pf_eig; + Eigen::Matrix baseline_work; + Eigen::Matrix baseline_other; + + // clang-format off + baseline_home << + 0.4413, 0.4504, 1.2383, 0.8033, 0.0494, 0.0017, + 0.0485, 0.7616, 0.6532, 1.1614, 0.0256, 0.0013, + 0.1800, 0.1795, 0.8806, 0.6413, 0.0429, 0.0032, + 0.0495, 0.2639, 0.5189, 0.8277, 0.0679, 0.0014, + 0.0087, 0.0394, 0.1417, 0.3834, 0.7064, 0.0447, + 0.0292, 0.0648, 0.1248, 0.4179, 0.3497, 0.1544; + + baseline_school_pf_eig << + 2.9964, 0.2501, 0.9132, 0.2509, 0.0000, 0.0000, + 0.2210, 1.9155, 0.2574, 0.2863, 0.0070, 0.0000, + 0.0324, 0.3598, 1.2613, 0.1854, 0.0041, 0.0000, + 0.1043, 0.4794, 1.1886, 0.1836, 0.0052, 0.0022, + 0.0000, 0.1150, 0.0000, 0.0000, 0.0168, 0.0000, + 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000; + + baseline_other << + 0.5170, 0.3997, 0.7957, 0.9958, 0.3239, 0.0428, + 0.0632, 0.9121, 0.3254, 0.4731, 0.2355, 0.0148, + 0.0336, 0.1604, 1.7529, 0.8622, 0.1440, 0.0077, + 0.0204, 0.1444, 0.5738, 1.2127, 0.3433, 0.0178, + 0.0371, 0.0393, 0.4171, 0.9666, 0.7495, 0.0257, + 0.0791, 0.0800, 0.3480, 0.5588, 0.2769, 0.0180; + + baseline_work << + 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, + 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, + 0.0000, 0.0127, 1.7570, 1.6050, 0.0133, 0.0000, + 0.0000, 0.0020, 1.0311, 2.3166, 0.0098, 0.0000, + 0.0000, 0.0002, 0.0194, 0.0325, 0.0003, 0.0000, + 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000; + // clang-format on + + Eigen::Matrix minimum_home; + Eigen::Matrix minimum_school_pf_eig; + Eigen::Matrix minimum_work; + Eigen::Matrix minimum_other; + + minimum_home.setZero(); + minimum_school_pf_eig.setZero(); + minimum_work.setZero(); + minimum_other.setZero(); + + mio::ContactMatrixGroup contact_matrices(num_contact_locations, num_age_groups); + + contact_matrices[0].get_baseline() = baseline_home; + contact_matrices[1].get_baseline() = baseline_school_pf_eig; + contact_matrices[2].get_baseline() = baseline_work; + contact_matrices[3].get_baseline() = baseline_other; + + contact_matrices[0].get_minimum() = minimum_home; + contact_matrices[1].get_minimum() = minimum_school_pf_eig; + contact_matrices[2].get_minimum() = minimum_work; + contact_matrices[3].get_minimum() = minimum_other; + + params.template get>() = + mio::UncertainContactMatrix(std::move(contact_matrices)); + + // ----------------------------- // + // Define the initial population // + // clang-format off + size_t num_infection_states = static_cast(mio::osecirvvs::InfectionState::Count); + Eigen::Matrix population(num_age_groups, num_infection_states); + population << + /* ------------- */ + /* AgeGroup: 0-4 */ + 2753030.497430, + 0.0, // Susceptible + 11.708820, 0.0, 2.212815, // Exposed + 7.578866, 0.0, 4.084619, 0.0, 0.0, 0.0, // Infected No Symptoms + 22.853384, 0.0, 2.146616, 0.0, 0.0, 0.0, // Infected Symptoms + 0.150378, 0.0, 0.001865, // Infected Severe + 0.150441, 0.0, 0.005732, // Infected Critical + 1051258.609035, // Recovered + 64.000000, 0.0, 0.0, // Dead + /* ------------- */ + /* AgeGroup: 5-14 */ + 3326384.779416, 0.0, // Susceptible + 5.015423, 0.0, 3.287046, // Exposed + 2.945401, 0.0, 5.907620, 0.0, 0.0, 0.0, // Infected No Symptoms + 5.168372, 0.0, 1.831628, 0.0, 0.0, 0.0, // Infected Symptoms + 0.023957, 0.0, 0.001460, // Infected Severe + 0.138450, 0.0, 0.017723, // Infected Critical + 4506996.883505, // Recovered + 45.000000, 0.0, 0.0, // Dead + /* ------------- */ + /* AgeGroup: 15-34 */ + 7474316.734280, 0.0, // Susceptible + 15.185091, 0.0, 9.963930, // Exposed + 10.119921, 0.0, 19.347768, 0.0, 0.0, 0.0, // Infected No Symptoms + 27.662806, 0.0, 9.908623, 0.0, 0.0, 0.0, // Infected Symptoms + 0.839570, 0.0, 0.044404, // Infected Severe + 0.349355, 0.0, 0.046282, // Infected Critical + 11247005.297970, // Recovered + 510.0, 0.0, 0.0, // Dead + /* ------------- */ + /* AgeGroup: 35-59 */ + 13079641.709818, 0.0, // Susceptible + 28.087697, 0.0, 14.523549, // Exposed + 23.028794, 0.0, 36.771572, 0.0, 0.0, 0.0, // Infected No Symptoms + 74.950401, 0.0, 21.335313, 0.0, 0.0, 0.0, // Infected Symptoms + 7.139552, 0.0, 0.287452, // Infected Severe + 2.329300, 0.0, 0.231930, // Infected Critical + 15095393.104624, // Recovered + 8936.000000, 0.0, 0.0, // Dead + /* ------------- */ + /* AgeGroup: 60-79 */ + 12370660.206757, 0.0, // Susceptible + 66.165784, 0.0, 13.006396, // Exposed + 51.098574, 0.0, 30.752489, 0.0, 0.0, 0.0, // Infected No Symptoms + 162.557513, 0.0, 16.728201, 0.0, 0.0, 0.0, // Infected Symptoms + 38.893189, 0.0, 0.564814, // Infected Severe + 13.248394, 0.0, 0.494789, // Infected Critical + 5239420.683100, // Recovered + 56248.571429, 0.0, 0.0, // Dead + /* ------------- */ + /* AgeGroup: 80+ */ + 5614932.848062, 0.0, // Susceptible + 61.140294, 0.0, 9.149288, // Exposed + 46.298860, 0.0, 21.109925, 0.0, 0.0, 0.0, // Infected No Symptoms + 137.630764, 0.0, 11.369236, 0.0, 0.0, 0.0, // Infected Symptoms + 59.433799, 0.0, 0.673454, // Infected Severe + 24.245028, 0.0, 0.742577, // Infected Critical + 1695860.958713, // Recovered + 121825.428571, 0.0, 0.0; // Dead + // clang-format on + + for (size_t column = 0; column < static_cast(population.cols()); column++) { + for (size_t row = 0; row < static_cast(population.rows()); row++) { + mio::Index index{mio::AgeGroup(row), + mio::osecirvvs::InfectionState(column)}; + model.populations[index] = population(row, column); + } + } + + model.apply_constraints(); + + // ---------------------- // + // Define the graph model // + mio::Graph>, mio::MobilityEdge> graph_simulation; + + graph_simulation.add_node(0, model); + + return graph_simulation; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_settings/optimization_settings.cpp b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_settings/optimization_settings.cpp new file mode 100644 index 0000000000..eba826d357 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_settings/optimization_settings.cpp @@ -0,0 +1,145 @@ +#include "optimization_settings.h" + +#include "../control_parameters/damping_controls.h" +#include "../constraints/update_constraints.h" + +SecirvvsOptimization::SecirvvsOptimization(const OptimizationModel& optimization_model, size_t num_control_intervals, + size_t pc_resolution, bool random_start, IntegratorType integrator_type, + size_t integrator_resolution, ADType ad_eval_f, ADType ad_eval_jac, + std::vector control_parameters, + std::vector path_constraints, + std::vector terminal_constraints, ControlActivation activation) + : m_optimization_model(optimization_model) + , m_t0(0.0) + , m_tmax(m_optimization_model.simulation_days()) + , m_num_control_intervals(num_control_intervals) + , m_pc_resolution(pc_resolution) + , m_random_start(random_start) + , m_integrator_type(integrator_type) + , m_integrator_resolution(integrator_resolution) + , m_ad_eval_f(ad_eval_f) + , m_ad_eval_jac(ad_eval_jac) + , m_control_parameters(control_parameters) + , m_path_constraints(path_constraints) + , m_terminal_constraints(terminal_constraints) + , m_activation_function(ControlActivationFunction(activation)) +{ + // Number of intervals to check for path constraints (each day). + m_num_intervals = m_num_control_intervals * m_pc_resolution; + // Time step size used in the integrator. + m_dt = (m_tmax - m_t0) / (m_num_intervals * m_integrator_resolution); + // Number of parameter and constraints. + m_num_control_parameters = control_parameters.size(); + m_num_path_constraints = path_constraints.size(); + m_num_terminal_constraints = terminal_constraints.size(); +} + +const OptimizationModel& SecirvvsOptimization::optimization_model() const +{ + return this->m_optimization_model; +} + +double SecirvvsOptimization::t0() const +{ + return m_t0; +} + +double SecirvvsOptimization::tmax() const +{ + return m_tmax; +} + +size_t SecirvvsOptimization::num_control_intervals() const +{ + return this->m_num_control_intervals; +} + +size_t SecirvvsOptimization::pc_resolution() const +{ + return this->m_pc_resolution; +} + +bool SecirvvsOptimization::random_start() const +{ + return this->m_random_start; +} + +IntegratorType SecirvvsOptimization::integrator_type() const +{ + return this->m_integrator_type; +} + +size_t SecirvvsOptimization::integrator_resolution() const +{ + return this->m_integrator_resolution; +} + +ADType SecirvvsOptimization::ad_eval_f() const +{ + return this->m_ad_eval_f; +} + +ADType SecirvvsOptimization::ad_eval_jac() const +{ + return this->m_ad_eval_jac; +} + +std::vector& SecirvvsOptimization::control_parameters() +{ + return this->m_control_parameters; +} + +const std::vector& SecirvvsOptimization::control_parameters() const +{ + return this->m_control_parameters; +} + +std::vector& SecirvvsOptimization::path_constraints() +{ + return this->m_path_constraints; +} + +const std::vector& SecirvvsOptimization::path_constraints() const +{ + return this->m_path_constraints; +} + +std::vector& SecirvvsOptimization::terminal_constraints() +{ + return this->m_terminal_constraints; +} + +const std::vector& SecirvvsOptimization::terminal_constraints() const +{ + return this->m_terminal_constraints; +} + +ControlActivationFunction SecirvvsOptimization::activation_function() const +{ + return this->m_activation_function; +} + +size_t SecirvvsOptimization::num_intervals() const +{ + return this->m_num_intervals; +} + +double SecirvvsOptimization::dt() const +{ + return this->m_dt; +} + +size_t SecirvvsOptimization::num_control_parameters() const +{ + return this->m_num_control_parameters; +} + +size_t SecirvvsOptimization::num_path_constraints() const +{ + return this->m_num_path_constraints; +} + +size_t SecirvvsOptimization::num_terminal_constraints() const +{ + return this->m_num_terminal_constraints; +} diff --git a/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_settings/optimization_settings.h b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_settings/optimization_settings.h new file mode 100644 index 0000000000..79b1134ce8 --- /dev/null +++ b/cpp/tools/optimal_control/secirvvs_ESID_dampings/optimization_settings/optimization_settings.h @@ -0,0 +1,80 @@ +#pragma once + +#include "../optimization_model/optimization_model.h" +#include "../control_parameters/control_parameters.h" +#include "../control_parameters/control_activation.h" +#include "../constraints/constraints.h" +#include "../helpers/integrator_selector.h" +#include "../helpers/ad_type.h" + +#include +#include +#include +#include +#include + +/* Stores the optimization model and the settings used in optimal control */ +class SecirvvsOptimization +{ +public: + SecirvvsOptimization(const OptimizationModel& optimization_model, size_t num_control_intervals, + size_t pc_resolution, bool random_start, IntegratorType integrator_type, + size_t integrator_resolution, ADType ad_eval_f, ADType ad_eval_jac, + std::vector control_parameters, + std::vector path_constraints, std::vector terminal_constraints, + ControlActivation activation); + + const OptimizationModel& optimization_model() const; + double t0() const; + double tmax() const; + + size_t num_control_intervals() const; // Number of control intervals + size_t pc_resolution() const; // Path constraint resolution + bool random_start() const; // Whether to randomize initial control values + IntegratorType integrator_type() const; // Type of integrator used + size_t integrator_resolution() const; // Resolution of the integrator + ADType ad_eval_f() const; // AD type for function evaluation + ADType ad_eval_jac() const; // AD type for Jacobian evaluation + + // Control parameters used in the optimization + std::vector& control_parameters(); + const std::vector& control_parameters() const; + // Path constraints used in the optimization + std::vector& path_constraints(); + const std::vector& path_constraints() const; + // Terminal constraints used in the optimization + std::vector& terminal_constraints(); + const std::vector& terminal_constraints() const; + + // Mapping f:[0,1]->[0,1] for control parameters + ControlActivationFunction activation_function() const; + + size_t num_intervals() const; // Number of intervals to check for path constraints (each day). + double dt() const; // Time step size used in the integrator. + size_t num_control_parameters() const; // Number of control parameters. + size_t num_path_constraints() const; // Number of path constraints. + size_t num_terminal_constraints() const; // Number of terminal constraints. + +private: + /* Constructor members */ + const OptimizationModel& m_optimization_model; + double m_t0; + double m_tmax; + size_t m_num_control_intervals; + size_t m_pc_resolution; + bool m_random_start; + IntegratorType m_integrator_type; + size_t m_integrator_resolution; + ADType m_ad_eval_f; + ADType m_ad_eval_jac; + std::vector m_control_parameters; + std::vector m_path_constraints; + std::vector m_terminal_constraints; + ControlActivationFunction m_activation_function; + /* Helper members */ + size_t m_num_intervals; + double m_dt; + size_t m_num_control_parameters; + size_t m_num_path_constraints; + size_t m_num_terminal_constraints; +};