From dc1ec15fc6fec59cc4c70fdb719bf731d7a57c9c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 23 Apr 2024 18:28:16 +0200
Subject: [PATCH] Refactor output for complementarity conditions
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The vector of lower and upper bounds is now provided by
{static,dynamic}_complementarity_conditions.m (which accepts parameters as an
argument, in preparation for the possibility of having parameters in
complementarity conditions).

The index that denotes the reordering of equation residuals (so that the
residual of an equation appears at the index of the endogenous appearing in the
complementarity condition) is now in “M_.{static,dynamic}_mcp_equations_reordering”.

Constraints declared through the ramsey_constraints block are also incorporated
in this new interface (and therefore “M_.ramsey_model_constraints” has been
removed).
---
 src/ComputingTasks.cc |  84 ------------------------------
 src/ComputingTasks.hh |  24 +--------
 src/DynamicModel.cc   |  81 ++++++++++++++++++++++-------
 src/DynamicModel.hh   |  22 +++++---
 src/DynareBison.yy    |  29 +++++------
 src/ExprNode.cc       |  39 +++++++++++++-
 src/ExprNode.hh       |   6 +++
 src/ModFile.cc        |  13 +++--
 src/ModFile.hh        |   7 ++-
 src/ModelTree.cc      |  86 ++++++++++++++++++++++++++++---
 src/ModelTree.hh      |  80 ++++++++++++++++++++++++++---
 src/ParsingDriver.cc  | 116 +++++++++++++++++++++++++-----------------
 src/ParsingDriver.hh  |  17 ++-----
 src/Statement.hh      |   2 -
 src/StaticModel.cc    |  14 ++++-
 15 files changed, 389 insertions(+), 231 deletions(-)

diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index f278bf6c..98b8d3a2 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -508,90 +508,6 @@ RamseyModelStatement::writeJsonOutput(ostream& output) const
   output << "}";
 }
 
-RamseyConstraintsStatement::RamseyConstraintsStatement(const SymbolTable& symbol_table_arg,
-                                                       constraints_t constraints_arg) :
-    symbol_table {symbol_table_arg}, constraints {move(constraints_arg)}
-{
-}
-
-void
-RamseyConstraintsStatement::checkPass(ModFileStructure& mod_file_struct,
-                                      [[maybe_unused]] WarningConsolidation& warnings)
-{
-  mod_file_struct.ramsey_constraints_present = true;
-}
-
-void
-RamseyConstraintsStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
-                                        [[maybe_unused]] bool minimal_workspace) const
-{
-  output << "M_.ramsey_model_constraints = {" << endl;
-  for (bool printed_something {false}; const auto& it : constraints)
-    {
-      if (exchange(printed_something, true))
-        output << ", ";
-      output << "{" << it.endo + 1 << ", '";
-      switch (it.code)
-        {
-        case BinaryOpcode::less:
-          output << '<';
-          break;
-        case BinaryOpcode::greater:
-          output << '>';
-          break;
-        case BinaryOpcode::lessEqual:
-          output << "<=";
-          break;
-        case BinaryOpcode::greaterEqual:
-          output << ">=";
-          break;
-        default:
-          cerr << "Ramsey constraints: this shouldn't happen." << endl;
-          exit(EXIT_FAILURE);
-        }
-      output << "', '";
-      it.expression->writeOutput(output);
-      output << "'}" << endl;
-    }
-  output << "};" << endl;
-}
-
-void
-RamseyConstraintsStatement::writeJsonOutput(ostream& output) const
-{
-  output << R"({"statementName": "ramsey_constraints")"
-         << R"(, "ramsey_model_constraints": [)" << endl;
-  for (bool printed_something {false}; const auto& it : constraints)
-    {
-      if (exchange(printed_something, true))
-        output << ", ";
-      output << R"({"constraint": ")" << symbol_table.getName(it.endo) << " ";
-      switch (it.code)
-        {
-        case BinaryOpcode::less:
-          output << '<';
-          break;
-        case BinaryOpcode::greater:
-          output << '>';
-          break;
-        case BinaryOpcode::lessEqual:
-          output << "<=";
-          break;
-        case BinaryOpcode::greaterEqual:
-          output << ">=";
-          break;
-        default:
-          cerr << "Ramsey constraints: this shouldn't happen." << endl;
-          exit(EXIT_FAILURE);
-        }
-      output << " ";
-      it.expression->writeJsonOutput(output, {}, {});
-      output << R"("})" << endl;
-    }
-  output << "]" << endl;
-  output << "}";
-}
-
 RamseyPolicyStatement::RamseyPolicyStatement(SymbolList symbol_list_arg,
                                              OptionsList options_list_arg,
                                              const SymbolTable& symbol_table_arg) :
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 8909a8c7..747bb0a5 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -179,28 +179,6 @@ public:
   void writeJsonOutput(ostream& output) const override;
 };
 
-class RamseyConstraintsStatement : public Statement
-{
-public:
-  struct Constraint
-  {
-    int endo;
-    BinaryOpcode code;
-    expr_t expression;
-  };
-  using constraints_t = vector<Constraint>;
-
-private:
-  const SymbolTable& symbol_table;
-  const constraints_t constraints;
-
-public:
-  RamseyConstraintsStatement(const SymbolTable& symbol_table_arg, constraints_t constraints_arg);
-  void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
-  void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
-  void writeJsonOutput(ostream& output) const override;
-};
-
 class RamseyPolicyStatement : public Statement
 {
 private:
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 50b02216..76323b77 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -34,10 +34,18 @@
 void
 DynamicModel::copyHelper(const DynamicModel& m)
 {
-  auto f = [this](const ExprNode* e) { return e->clone(*this); };
+  auto f = [this](const ExprNode* e) { return e ? e->clone(*this) : nullptr; };
 
   for (const auto& it : m.static_only_equations)
     static_only_equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
+  for (const auto& it : m.static_only_complementarity_conditions)
+    if (it)
+      {
+        const auto& [symb_id, lb, ub] = *it;
+        static_only_complementarity_conditions.emplace_back(in_place, symb_id, f(lb), f(ub));
+      }
+    else
+      static_only_complementarity_conditions.emplace_back(nullopt);
 }
 
 DynamicModel::DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
@@ -104,6 +112,7 @@ DynamicModel::operator=(const DynamicModel& m)
 
   static_only_equations_lineno = m.static_only_equations_lineno;
   static_only_equations_equation_tags = m.static_only_equations_equation_tags;
+  static_only_complementarity_conditions.clear();
   deriv_id_table = m.deriv_id_table;
   inv_deriv_id_table = m.inv_deriv_id_table;
   dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
@@ -651,11 +660,11 @@ DynamicModel::parseIncludeExcludeEquations(const string& inc_exc_option_value, b
 }
 
 vector<int>
-DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs,
-                                    bool excluded_vars_change_type,
-                                    vector<BinaryOpNode*>& all_equations,
-                                    vector<optional<int>>& all_equations_lineno,
-                                    EquationTags& all_equation_tags, bool static_equations) const
+DynamicModel::removeEquationsHelper(
+    set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
+    vector<BinaryOpNode*>& all_equations, vector<optional<int>>& all_equations_lineno,
+    vector<optional<tuple<int, expr_t, expr_t>>>& all_complementarity_conditions,
+    EquationTags& all_equation_tags, bool static_equations) const
 {
   if (all_equations.empty())
     return {};
@@ -686,6 +695,7 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
   // remove from equations, equations_lineno, equation_tags
   vector<BinaryOpNode*> new_equations;
   vector<optional<int>> new_equations_lineno;
+  vector<optional<tuple<int, expr_t, expr_t>>> new_complementarity_conditions;
   map<int, int> old_eqn_num_2_new;
   vector<int> excluded_vars;
   for (size_t i = 0; i < all_equations.size(); i++)
@@ -717,11 +727,13 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
         new_equations.emplace_back(all_equations[i]);
         old_eqn_num_2_new[i] = new_equations.size() - 1;
         new_equations_lineno.emplace_back(all_equations_lineno[i]);
+        new_complementarity_conditions.emplace_back(all_complementarity_conditions[i]);
       }
   int n_excl = all_equations.size() - new_equations.size();
 
   all_equations = new_equations;
   all_equations_lineno = new_equations_lineno;
+  all_complementarity_conditions = new_complementarity_conditions;
 
   all_equation_tags.erase(eqs_to_delete_by_number, old_eqn_num_2_new);
 
@@ -754,12 +766,13 @@ DynamicModel::removeEquations(const vector<map<string, string>>& listed_eqs_by_t
 
   vector<int> excluded_vars
       = removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type, equations,
-                              equations_lineno, equation_tags, false);
+                              equations_lineno, complementarity_conditions, equation_tags, false);
 
   // Ignore output because variables are not excluded when equations marked 'static' are excluded
   removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type,
                         static_only_equations, static_only_equations_lineno,
-                        static_only_equations_equation_tags, true);
+                        static_only_complementarity_conditions, static_only_equations_equation_tags,
+                        true);
 
   if (!listed_eqs_by_tag2.empty())
     {
@@ -1134,6 +1147,11 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
       output << "'; " << endl;
     }
   output << "};" << endl;
+
+  output << "M_.dynamic_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
 }
 
 void
@@ -2507,6 +2525,8 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder,
       cerr << "ERROR: Block decomposition requested but failed." << endl;
       exit(EXIT_FAILURE);
     }
+
+  computeMCPEquationsReordering();
 }
 
 void
@@ -2810,6 +2830,8 @@ DynamicModel::writeDynamicFile(const string& basename, bool use_dll, const strin
 
   writeSetAuxiliaryVariablesFile<true>(basename, julia);
 
+  writeComplementarityConditionsFile<true>(basename);
+
   // Support for model debugging
   if (!julia)
     writeDebugModelMFiles<true>(basename);
@@ -2821,6 +2843,7 @@ DynamicModel::clearEquations()
   equations.clear();
   equations_lineno.clear();
   equation_tags.clear();
+  complementarity_conditions.clear();
 }
 
 void
@@ -2828,14 +2851,26 @@ DynamicModel::replaceMyEquations(DynamicModel& dynamic_model) const
 {
   dynamic_model.clearEquations();
 
+  auto clone_if_not_null = [&](expr_t e) { return e ? e->clone(dynamic_model) : nullptr; };
+
   for (size_t i = 0; i < equations.size(); i++)
-    dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i]);
+    {
+      optional<tuple<int, expr_t, expr_t>> cc_clone;
+      if (complementarity_conditions[i])
+        {
+          auto& [symb_id, lower_bound, upper_bound] = *complementarity_conditions[i];
+          cc_clone = {symb_id, clone_if_not_null(lower_bound), clone_if_not_null(upper_bound)};
+        }
+      dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i],
+                                move(cc_clone));
+    }
 
   dynamic_model.equation_tags = equation_tags;
 }
 
 int
-DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
+DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model,
+                                      map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints)
 {
   cout << "Ramsey Problem: added " << equations.size() << " multipliers." << endl;
 
@@ -2892,9 +2927,10 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
                       lagrangian);
       }
 
-  // Save line numbers and tags, see below
+  // Save line numbers, tags and complementarity conditions, see below
   auto old_equations_lineno = equations_lineno;
   auto old_equation_tags = equation_tags;
+  auto old_complementarity_conditions = complementarity_conditions;
 
   // Prepare derivation of the Lagrangian
   clearEquations();
@@ -2911,6 +2947,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
   vector<expr_t> neweqs;
   vector<optional<int>> neweqs_lineno;
   map<int, map<string, string>> neweqs_tags;
+  map<int, optional<tuple<int, expr_t, expr_t>>> new_complementarity_conditions;
   int orig_endo_nbr {0};
   for (auto& [symb_id_and_lag, deriv_id] : deriv_id_table)
     {
@@ -2924,11 +2961,19 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
               // This is a derivative w.r.t. a Lagrange multiplier
               neweqs_lineno.push_back(old_equations_lineno[*i]);
               neweqs_tags[neweqs.size() - 1] = old_equation_tags.getTagsByEqn(*i);
+              new_complementarity_conditions.emplace(neweqs.size() - 1,
+                                                     old_complementarity_conditions.at(*i));
             }
           else
             {
               orig_endo_nbr++;
               neweqs_lineno.emplace_back(nullopt);
+              if (cloned_ramsey_constraints.contains(symb_id))
+                {
+                  auto& [lower_bound, upper_bound] = cloned_ramsey_constraints.at(symb_id);
+                  new_complementarity_conditions.emplace(neweqs.size() - 1,
+                                                         tuple {symb_id, lower_bound, upper_bound});
+                }
             }
         }
     }
@@ -2936,7 +2981,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
   // Overwrite equations with the Lagrangian derivatives
   clearEquations();
   for (size_t i = 0; i < neweqs.size(); i++)
-    addEquation(neweqs[i], neweqs_lineno[i], neweqs_tags[i]);
+    addEquation(neweqs[i], neweqs_lineno[i], new_complementarity_conditions[i], neweqs_tags[i]);
 
   return orig_endo_nbr;
 }
@@ -3773,6 +3818,7 @@ DynamicModel::fillEvalContext(eval_context_t& eval_context) const
 
 void
 DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
+                                    optional<tuple<int, expr_t, expr_t>> complementarity_condition,
                                     map<string, string> eq_tags)
 {
   auto beq = dynamic_cast<BinaryOpNode*>(eq);
@@ -3781,6 +3827,7 @@ DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
   static_only_equations_equation_tags.add(static_only_equations.size(), move(eq_tags));
   static_only_equations.push_back(beq);
   static_only_equations_lineno.push_back(lineno);
+  static_only_complementarity_conditions.push_back(move(complementarity_condition));
 }
 
 size_t
@@ -3834,7 +3881,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
     {
       auto eq_tags_dynamic = eq_tags;
       eq_tags_dynamic["dynamic"] = "";
-      addEquation(AddEqual(term, Zero), lineno, eq_tags_dynamic);
+      addEquation(AddEqual(term, Zero), lineno, nullopt, eq_tags_dynamic);
     }
 
   // Create or update the static equation (corresponding to the pure relax regime)
@@ -3856,7 +3903,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
       else
         {
           eq_tags["static"] = "";
-          addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, move(eq_tags));
+          addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, nullopt, move(eq_tags));
         }
     }
 }
@@ -4171,8 +4218,8 @@ DynamicModel::simplifyEquations()
 {
   size_t last_subst_table_size = 0;
   map<VariableNode*, NumConstNode*> subst_table;
-  // Equations with “mcp” tag are excluded, see dynare#1697
-  findConstantEquationsWithoutMcpTag(subst_table);
+  // Equations with a complementarity condition are excluded, see dynare#1697
+  findConstantEquationsWithoutComplementarityCondition(subst_table);
   while (subst_table.size() != last_subst_table_size)
     {
       last_subst_table_size = subst_table.size();
@@ -4183,7 +4230,7 @@ DynamicModel::simplifyEquations()
       for (auto& equation : static_only_equations)
         equation = dynamic_cast<BinaryOpNode*>(equation->replaceVarsInEquation(subst_table));
       subst_table.clear();
-      findConstantEquationsWithoutMcpTag(subst_table);
+      findConstantEquationsWithoutComplementarityCondition(subst_table);
     }
 }
 
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 262d7599..a07f4b6f 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -92,6 +92,9 @@ private:
   //! Stores the equation tags of equations declared as [static]
   EquationTags static_only_equations_equation_tags;
 
+  // Complementarity conditions of equations declared as [static]
+  vector<optional<tuple<int, expr_t, expr_t>>> static_only_complementarity_conditions;
+
   using deriv_id_table_t = map<pair<int, int>, int>;
   //! Maps a pair (symbol_id, lag) to a deriv ID
   deriv_id_table_t deriv_id_table;
@@ -279,11 +282,11 @@ private:
 
      Returns a list of excluded variables (empty if
      excluded_vars_change_type=false) */
-  vector<int> removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs,
-                                    bool excluded_vars_change_type,
-                                    vector<BinaryOpNode*>& all_equations,
-                                    vector<optional<int>>& all_equations_lineno,
-                                    EquationTags& all_equation_tags, bool static_equations) const;
+  vector<int> removeEquationsHelper(
+      set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
+      vector<BinaryOpNode*>& all_equations, vector<optional<int>>& all_equations_lineno,
+      vector<optional<tuple<int, expr_t, expr_t>>>& all_complementarity_conditions,
+      EquationTags& all_equation_tags, bool static_equations) const;
 
   //! Compute autoregressive matrices of trend component models
   /* The algorithm uses matching rules over expression trees. It cannot handle
@@ -455,7 +458,8 @@ public:
      Returns the number of optimality FOCs, which is by construction equal to
      the number of endogenous before adding the Lagrange multipliers
      (internally called ramsey_endo_nbr). */
-  int computeRamseyPolicyFOCs(const StaticModel& static_model);
+  int computeRamseyPolicyFOCs(const StaticModel& _model,
+                              map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints);
 
   //! Clears all equations
   void clearEquations();
@@ -464,7 +468,9 @@ public:
   void replaceMyEquations(DynamicModel& dynamic_model) const;
 
   //! Adds an equation marked as [static]
-  void addStaticOnlyEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags);
+  void addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
+                             optional<tuple<int, expr_t, expr_t>> complementarity_condition,
+                             map<string, string> eq_tags);
 
   //! Returns number of static only equations
   size_t staticOnlyEquationsNbr() const;
@@ -708,7 +714,7 @@ public:
   getStaticOnlyEquationsInfo() const
   {
     return tuple {static_only_equations, static_only_equations_lineno,
-                  static_only_equations_equation_tags};
+                  static_only_complementarity_conditions, static_only_equations_equation_tags};
   };
 
   //! Returns true if a parameter was used in the model block with a lead or lag
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 4aed10f5..b014b707 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -242,7 +242,7 @@ str_tolower(string s)
 %type <map<string, string>> tag_pair_list
 %type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt
 %type <vector<pair<int, int>>> period_list
-%type <vector<expr_t>> matched_moments_list value_list
+%type <vector<expr_t>> matched_moments_list value_list ramsey_constraints_list
 %type <tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> occbin_constraints_regime
 %type <vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>>> occbin_constraints_regimes_list
 %type <map<string, expr_t>> occbin_constraints_regime_options_list
@@ -2616,23 +2616,20 @@ ramsey_policy : RAMSEY_POLICY ';'
                 { driver.ramsey_policy($5); }
               ;
 
-ramsey_constraints : RAMSEY_CONSTRAINTS ';' ramsey_constraints_list END ';'
-                     { driver.add_ramsey_constraints_statement(); }
+ramsey_constraints : RAMSEY_CONSTRAINTS ';'
+                     { driver.begin_ramsey_constraints(); }
+                     ramsey_constraints_list END ';'
+                     { driver.end_ramsey_constraints($4); }
 		   ;
 
-ramsey_constraints_list : ramsey_constraints_list ramsey_constraint
-                 | ramsey_constraint
-		 ;
-
-ramsey_constraint : NAME  LESS expression ';'
-                    { driver.ramsey_constraint_add_less($1, $3); }
-		  | NAME  GREATER  expression ';'
-                    { driver.ramsey_constraint_add_greater($1, $3); }
-		  | NAME  LESS_EQUAL expression ';'
-                    { driver.ramsey_constraint_add_less_equal($1, $3); }
-		  | NAME  GREATER_EQUAL  expression ';'
-                    { driver.ramsey_constraint_add_greater_equal($1, $3); }
-		  ;
+ramsey_constraints_list : ramsey_constraints_list hand_side ';'
+                          {
+                            $$ = $1;
+                            $$.push_back($2);
+                          }
+                        | hand_side ';'
+                          { $$ = { $1 }; }
+                        ;
 
 evaluate_planner_objective : EVALUATE_PLANNER_OBJECTIVE ';'
                              { driver.evaluate_planner_objective(); }
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index b3b31d71..44923c17 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -9370,3 +9370,40 @@ ExprNode::toString() const
   writeJsonOutput(ss, {}, {});
   return ss.str();
 }
+
+tuple<int, expr_t, expr_t>
+ExprNode::matchComplementarityCondition() const
+{
+  throw MatchFailureException {"This expression is not an inequality"};
+}
+
+tuple<int, expr_t, expr_t>
+BinaryOpNode::matchComplementarityCondition() const
+{
+  bool is_lower_bound {[&] {
+    switch (op_code)
+      {
+      case BinaryOpcode::less:
+      case BinaryOpcode::lessEqual:
+        return false;
+      case BinaryOpcode::greater:
+      case BinaryOpcode::greaterEqual:
+        return true;
+      default:
+        throw MatchFailureException {"This expression is not an inequality"};
+      }
+  }()};
+
+  auto* varg = dynamic_cast<VariableNode*>(arg1);
+  if (!varg)
+    throw MatchFailureException {"Left-hand side is not a variable"};
+  if (varg->lag != 0)
+    throw MatchFailureException {"Left-hand side variable must not have a lead or a lag"};
+  if (datatree.symbol_table.getType(varg->symb_id) != SymbolType::endogenous)
+    throw MatchFailureException {"Left-hand side is not an endogenous variable"};
+
+  if (!arg2->isConstant())
+    throw MatchFailureException {"Right-hand side is not a constant"};
+
+  return {varg->symb_id, is_lower_bound ? arg2 : nullptr, is_lower_bound ? nullptr : arg2};
+}
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index a958a22d..c4d73c28 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -942,6 +942,11 @@ public:
 
   // Substitutes orig_symb_id(±l) with exp(aux_symb_id(±l)) (used for “var(log)”)
   [[nodiscard]] virtual expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const = 0;
+
+  /* Matches an expression that constitutes a complementarity condition.
+     If successful, returns a triplet (endo_symb_id, lower_bound, upper_bound).
+     Otherwise, throws a MatchFailureException. */
+  [[nodiscard]] virtual tuple<int, expr_t, expr_t> matchComplementarityCondition() const;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -1499,6 +1504,7 @@ public:
                           vector<int>& powers) const override;
   [[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
   [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
 };
 
 //! Trinary operator node
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 6ab8b267..3c01dbd1 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2006-2023 Dynare Team
+ * Copyright © 2006-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -167,7 +167,7 @@ ModFile::checkPass(bool nostrict, bool stochastic)
       exit(EXIT_FAILURE);
     }
 
-  if (mod_file_struct.ramsey_constraints_present && !mod_file_struct.ramsey_model_present)
+  if (!ramsey_constraints.empty() && !mod_file_struct.ramsey_model_present)
     {
       cerr << "ERROR: A ramsey_constraints block requires the presence of a ramsey_model or "
               "ramsey_policy statement"
@@ -535,8 +535,15 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
           symbol_table, num_constants, external_functions_table, trend_component_model_table,
           var_model_table};
       ramsey_FOC_equations_dynamic_model = dynamic_model;
+      auto clone_if_not_null
+          = [&](expr_t e) { return e ? e->clone(ramsey_FOC_equations_dynamic_model) : nullptr; };
+      map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints;
+      for (const auto& [symb_id, bounds] : ramsey_constraints)
+        cloned_ramsey_constraints.try_emplace(symb_id, clone_if_not_null(bounds.first),
+                                              clone_if_not_null(bounds.second));
       mod_file_struct.ramsey_orig_endo_nbr
-          = ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective);
+          = ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective,
+                                                                       cloned_ramsey_constraints);
       ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
     }
 
diff --git a/src/ModFile.hh b/src/ModFile.hh
index 440fc01f..22d50f9a 100644
--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2006-2023 Dynare Team
+ * Copyright © 2006-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -117,6 +117,11 @@ public:
   /*! (i.e. option parallel_local_files of model block) */
   vector<string> parallel_local_files;
 
+  /* Contents of the ramsey_constraints block.
+     Maps symb_id → (lower_bound, upper_bound).
+     NB: The two expr_t live in dynamic_model */
+  map<int, pair<expr_t, expr_t>> ramsey_constraints;
+
 private:
   //! List of statements
   vector<unique_ptr<Statement>> statements;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 41bed024..f2313f6f 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -34,6 +34,7 @@
 #endif
 
 #include <algorithm>
+#include <numeric>
 #include <regex>
 #include <utility>
 
@@ -51,13 +52,21 @@ vector<jthread> ModelTree::mex_compilation_workers;
 void
 ModelTree::copyHelper(const ModelTree& m)
 {
-  auto f = [this](expr_t e) { return e->clone(*this); };
+  auto f = [this](expr_t e) { return e ? e->clone(*this) : nullptr; };
 
   // Equations
   for (const auto& it : m.equations)
     equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
   for (const auto& it : m.aux_equations)
     aux_equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
+  for (const auto& it : m.complementarity_conditions)
+    if (it)
+      {
+        const auto& [symb_id, lb, ub] = *it;
+        complementarity_conditions.emplace_back(in_place, symb_id, f(lb), f(ub));
+      }
+    else
+      complementarity_conditions.emplace_back(nullopt);
 
   auto convert_deriv_map = [f](const map<vector<int>, expr_t>& dm) {
     map<vector<int>, expr_t> dm2;
@@ -184,6 +193,8 @@ ModelTree::operator=(const ModelTree& m)
   equations_lineno = m.equations_lineno;
   aux_equations.clear();
   equation_tags = m.equation_tags;
+  complementarity_conditions.clear();
+
   computed_derivs_order = m.computed_derivs_order;
   NNZDerivatives = m.NNZDerivatives;
 
@@ -1394,28 +1405,33 @@ ModelTree::writeLatexModelFile(const string& mod_basename, const string& latex_b
 }
 
 void
-ModelTree::addEquation(expr_t eq, const optional<int>& lineno)
+ModelTree::addEquation(expr_t eq, const optional<int>& lineno,
+                       optional<tuple<int, expr_t, expr_t>> complementarity_condition)
 {
   auto beq = dynamic_cast<BinaryOpNode*>(eq);
   assert(beq && beq->op_code == BinaryOpcode::equal);
 
   equations.push_back(beq);
   equations_lineno.push_back(lineno);
+  complementarity_conditions.push_back(move(complementarity_condition));
 }
 
 void
-ModelTree::findConstantEquationsWithoutMcpTag(map<VariableNode*, NumConstNode*>& subst_table) const
+ModelTree::findConstantEquationsWithoutComplementarityCondition(
+    map<VariableNode*, NumConstNode*>& subst_table) const
 {
   for (size_t i = 0; i < equations.size(); i++)
-    if (!equation_tags.exists(i, "mcp"))
+    if (!complementarity_conditions[i])
       equations[i]->findConstantEquations(subst_table);
 }
 
 void
-ModelTree::addEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags)
+ModelTree::addEquation(expr_t eq, const optional<int>& lineno,
+                       optional<tuple<int, expr_t, expr_t>> complementarity_condition,
+                       map<string, string> eq_tags)
 {
   equation_tags.add(equations.size(), move(eq_tags));
-  addEquation(eq, lineno);
+  addEquation(eq, lineno, move(complementarity_condition));
 }
 
 void
@@ -1593,6 +1609,27 @@ ModelTree::writeJsonModelEquations(ostream& output, bool residuals) const
               eqtags.clear();
             }
         }
+
+      if (complementarity_conditions[eq])
+        {
+          auto& [symb_id, lower_bound, upper_bound] = *complementarity_conditions[eq];
+          output << R"(, "complementarity_condition": {"variable": ")"
+                 << symbol_table.getName(symb_id) << '"';
+          if (lower_bound)
+            {
+              output << R"(, "lower_bound": ")";
+              lower_bound->writeJsonOutput(output, {}, {});
+              output << '"';
+            }
+          if (upper_bound)
+            {
+              output << R"(, "upper_bound": ")";
+              upper_bound->writeJsonOutput(output, {}, {});
+              output << '"';
+            }
+          output << "}";
+        }
+
       output << "}" << endl;
     }
   output << endl << "]" << endl;
@@ -2078,3 +2115,38 @@ ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType o
       output << ";" << endl;
     }
 }
+
+void
+ModelTree::computeMCPEquationsReordering()
+{
+  /* Optimal policy models (discretionary, or Ramsey before computing FOCs) do not have as many
+     equations as variables. Do not even try to compute the reordering. */
+  if (static_cast<int>(equations.size()) != symbol_table.endo_nbr())
+    return;
+
+  assert(equations.size() == complementarity_conditions.size());
+
+  mcp_equations_reordering.resize(equations.size());
+  iota(mcp_equations_reordering.begin(), mcp_equations_reordering.end(), 0);
+
+  set<int> endos;
+
+  for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
+    if (complementarity_conditions.at(eq))
+      {
+        int symb_id {get<0>(*complementarity_conditions[eq])};
+        auto [ignore, inserted] = endos.insert(symb_id);
+        if (!inserted)
+          {
+            cerr << "ERROR: variable " << symbol_table.getName(symb_id)
+                 << " appears in two complementarity conditions" << endl;
+            exit(EXIT_FAILURE);
+          }
+
+        int endo_id {symbol_table.getTypeSpecificID(symb_id)};
+
+        auto it = ranges::find(mcp_equations_reordering, eq);
+        assert(it != mcp_equations_reordering.end());
+        swap(mcp_equations_reordering[endo_id], *it);
+      }
+}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index d241ddff..512ec862 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -91,6 +91,9 @@ protected:
   vector<optional<int>> equations_lineno;
   //! Stores equation tags
   EquationTags equation_tags;
+  /* The tuple is: endogenous symbol ID, lower bound (possibly nullptr), upper bound (possibly
+     nullptr). */
+  vector<optional<tuple<int, expr_t, expr_t>>> complementarity_conditions;
   /*
    * ************** END **************
    */
@@ -271,6 +274,12 @@ protected:
      Same remark as above regarding blocks of type “evaluate”. */
   vector<vector<int>> blocks_jacobian_sparse_colptr;
 
+  /* Indices of reordered equations for use with an MCP solver. Contains a permutation of the
+     equation indices, so that reordered equations appear at the (type-specific) index of the
+     endogenous to which they are associated through the complementarity condition.
+     Also see computeMCPEquationsReordering() method. */
+  vector<int> mcp_equations_reordering;
+
   //! Computes derivatives
   /*! \param order the derivation order
       \param vars the derivation IDs w.r.t. which compute the derivatives */
@@ -284,6 +293,11 @@ protected:
 
   //! Computes temporary terms for the file containing parameters derivatives
   void computeParamsDerivativesTemporaryTerms();
+
+  /* Computes the mcp_equations_reordering vector.
+     Also checks that a variable does not appear as constrained in two different equations. */
+  void computeMCPEquationsReordering();
+
   //! Writes temporary terms
   template<ExprNodeOutputType output_type>
   void writeTemporaryTerms(const temporary_terms_t& tt, temporary_terms_t& temp_term_union,
@@ -424,6 +438,9 @@ protected:
   template<bool dynamic>
   void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const;
 
+  template<bool dynamic>
+  void writeComplementarityConditionsFile(const string& basename) const;
+
 private:
   //! Sparse matrix of double to store the values of the static Jacobian
   /*! First index is equation number, second index is endogenous type specific ID */
@@ -651,10 +668,14 @@ protected:
 public:
   //! Absolute value under which a number is considered to be zero
   double cutoff {1e-15};
-  //! Declare a node as an equation of the model; also give its line number
-  void addEquation(expr_t eq, const optional<int>& lineno);
+  /* Declare a node as an equation of the model; also give its line number and complementarity
+     condition */
+  void addEquation(expr_t eq, const optional<int>& lineno,
+                   optional<tuple<int, expr_t, expr_t>> complementarity_condition = nullopt);
   //! Declare a node as an equation of the model, also giving its tags
-  void addEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags);
+  void addEquation(expr_t eq, const optional<int>& lineno,
+                   optional<tuple<int, expr_t, expr_t>> complementarity_condition,
+                   map<string, string> eq_tags);
   //! Declare a node as an auxiliary equation of the model, adding it at the end of the list of
   //! auxiliary equations
   void addAuxEquation(expr_t eq);
@@ -671,9 +692,10 @@ public:
   /*! Reorder auxiliary variables so that they appear in recursive order in
       set_auxiliary_variables.m and dynamic_set_auxiliary_series.m */
   void reorderAuxiliaryEquations();
-  //! Find equations of the form “variable=constant”, excluding equations with “mcp” tag (see
-  //! dynare#1697)
-  void findConstantEquationsWithoutMcpTag(map<VariableNode*, NumConstNode*>& subst_table) const;
+  /* Find equations of the form “variable=constant”, excluding equations with a complementarity
+     condition (see dynare#1697) */
+  void findConstantEquationsWithoutComplementarityCondition(
+      map<VariableNode*, NumConstNode*>& subst_table) const;
   /* Given an expression, searches for the first equation that has exactly this
      expression on the LHS, and returns the RHS of that equation.
      If no such equation can be found, throws an ExprNode::MatchFailureExpression */
@@ -3133,4 +3155,50 @@ ModelTree::writeSetAuxiliaryVariablesFile(const string& basename, bool julia) co
     }
 }
 
+template<bool dynamic>
+void
+ModelTree::writeComplementarityConditionsFile(const string& basename) const
+{
+  // TODO: when C++20 support is complete, mark the following string constexpr
+  const string funcname {(dynamic ? "dynamic"s : "static"s) + "_complementarity_conditions"};
+  const filesystem::path filename {packageDir(basename) / (funcname + ".m")};
+  /* Can’t use matlabOutsideModel for output type, since it uses M_.
+     Static is ok even for the dynamic model, since there are no leads/lags. */
+  constexpr ExprNodeOutputType output_type {ExprNodeOutputType::matlabStaticModel};
+
+  ofstream output {filename, ios::out | ios::binary};
+  if (!output.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  output << "function [lb, ub] = " << funcname << "(params)" << endl
+         << "ub = inf(" << equations.size() << ",1);" << endl
+         << "lb = -ub;" << endl;
+
+  for (const auto& it : complementarity_conditions)
+    if (it)
+      {
+        const auto& [symb_id, lb, ub] = *it;
+        int endo_id {symbol_table.getTypeSpecificID(symb_id)};
+        if (lb)
+          {
+            output << "lb(" << endo_id + 1 << ")=";
+            lb->writeOutput(output, output_type);
+            output << ";" << endl;
+          }
+        if (ub)
+          {
+            output << "ub(" << endo_id + 1 << ")=";
+            ub->writeOutput(output, output_type);
+            output << ";" << endl;
+          }
+      }
+
+  output << "end" << endl;
+
+  output.close();
+}
+
 #endif
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 7f351fda..953535a7 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2635,6 +2635,9 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
     if (key == "endogenous")
       declare_or_change_type(SymbolType::endogenous, value);
 
+  optional<tuple<int, expr_t, expr_t>> matched_complementarity_condition;
+
+  // Handle legacy “mcp” tags, for backward compatibility
   if (eq_tags.contains("mcp"))
     {
       if (complementarity_condition)
@@ -2643,20 +2646,54 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
       else
         warning("Specifying complementarity conditions with the 'mcp' tag is obsolete. Please "
                 "consider switching to the new syntax using the perpendicular symbol.");
-    }
 
-  if (complementarity_condition)
-    {
-      if (auto bcomp = dynamic_cast<BinaryOpNode*>(complementarity_condition);
-          !(bcomp
-            && (bcomp->op_code == BinaryOpcode::less || bcomp->op_code == BinaryOpcode::lessEqual
-                || bcomp->op_code == BinaryOpcode::greater
-                || bcomp->op_code == BinaryOpcode::greaterEqual)))
-        error("The complementarity constraint must be an inequality.");
+      auto [var_name, is_lower_bound, constant] {[&]() -> tuple<string, bool, string> {
+        auto tagval = eq_tags.at("mcp");
+        if (auto less_pos = tagval.find('<'); less_pos != string::npos)
+          return {tagval.substr(0, less_pos), false, tagval.substr(less_pos + 1)};
+        else if (auto greater_pos = tagval.find('>'); greater_pos != string::npos)
+          return {tagval.substr(0, greater_pos), true, tagval.substr(greater_pos + 1)};
+        else
+          error("'mcp' tag does not contain an inequality");
+      }()};
+
+      int symb_id {[&] {
+        try
+          {
+            return mod_file->symbol_table.getID(var_name);
+          }
+        catch (SymbolTable::UnknownSymbolNameException&)
+          {
+            error("Left hand-side of expression in 'mcp' tag is not a variable");
+          }
+      }()};
+
+      if (mod_file->symbol_table.getType(symb_id) != SymbolType::endogenous)
+        error("Left hand-side of expression in 'mcp' tag is not an endogenous variable");
+
+      expr_t matched_constant {[&] {
+        char* str_end;
+        double d = strtod(constant.c_str(), &str_end);
+        if (str_end == constant.c_str())
+          error("Right hand-side of expression in 'mcp' tag should be a constant");
+        return data_tree->AddPossiblyNegativeConstant(d);
+      }()};
 
-      eq_tags.emplace("mcp", complementarity_condition->toString());
+      matched_complementarity_condition = {symb_id, is_lower_bound ? matched_constant : nullptr,
+                                           is_lower_bound ? nullptr : matched_constant};
     }
 
+  if (complementarity_condition)
+    try
+      {
+        matched_complementarity_condition
+            = complementarity_condition->matchComplementarityCondition();
+      }
+    catch (ExprNode::MatchFailureException& e)
+      {
+        error("Complementarity condition has an incorrect form: " + e.message);
+      }
+
   if (eq_tags.contains("static"))
     {
       // If the equation is tagged [static]
@@ -2664,7 +2701,8 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
         error("An equation tagged [static] cannot contain leads, lags, expectations or "
               "STEADY_STATE operators");
 
-      dynamic_model->addStaticOnlyEquation(id, location.begin.line, eq_tags);
+      dynamic_model->addStaticOnlyEquation(id, location.begin.line,
+                                           move(matched_complementarity_condition), eq_tags);
     }
   else if (eq_tags.contains("bind") || eq_tags.contains("relax"))
     {
@@ -2740,7 +2778,8 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
         }
     }
   else // General case
-    model_tree->addEquation(id, location.begin.line, move(eq_tags));
+    model_tree->addEquation(id, location.begin.line, move(matched_complementarity_condition),
+                            move(eq_tags));
 
   return id;
 }
@@ -3629,48 +3668,31 @@ ParsingDriver::prior_posterior_function(bool prior_func)
 }
 
 void
-ParsingDriver::add_ramsey_constraints_statement()
-{
-  mod_file->addStatement(
-      make_unique<RamseyConstraintsStatement>(mod_file->symbol_table, move(ramsey_constraints)));
-  ramsey_constraints.clear();
-}
-
-void
-ParsingDriver::ramsey_constraint_add_less(const string& name, const expr_t rhs)
-{
-  add_ramsey_constraint(name, BinaryOpcode::less, rhs);
-}
-
-void
-ParsingDriver::ramsey_constraint_add_greater(const string& name, const expr_t rhs)
+ParsingDriver::begin_ramsey_constraints()
 {
-  add_ramsey_constraint(name, BinaryOpcode::greater, rhs);
-}
-
-void
-ParsingDriver::ramsey_constraint_add_less_equal(const string& name, const expr_t rhs)
-{
-  add_ramsey_constraint(name, BinaryOpcode::lessEqual, rhs);
+  set_current_data_tree(&mod_file->dynamic_model);
 }
 
 void
-ParsingDriver::ramsey_constraint_add_greater_equal(const string& name, const expr_t rhs)
+ParsingDriver::end_ramsey_constraints(const vector<expr_t>& constraints)
 {
-  add_ramsey_constraint(name, BinaryOpcode::greaterEqual, rhs);
-}
+  for (expr_t c : constraints)
+    try
+      {
+        auto [symb_id, lower_bound, upper_bound] = c->matchComplementarityCondition();
 
-void
-ParsingDriver::add_ramsey_constraint(const string& name, BinaryOpcode op_code, const expr_t rhs)
-{
-  check_symbol_is_endogenous(name);
-  int symb_id = mod_file->symbol_table.getID(name);
+        auto [it, success]
+            = mod_file->ramsey_constraints.try_emplace(symb_id, lower_bound, upper_bound);
+        if (!success)
+          error("The ramsey_constraints block contains two constraints for variable "
+                + mod_file->symbol_table.getName(symb_id));
+      }
+    catch (ExprNode::MatchFailureException& e)
+      {
+        error("Ramsey constraint has an incorrect form: " + e.message);
+      }
 
-  RamseyConstraintsStatement::Constraint C;
-  C.endo = symb_id;
-  C.code = op_code;
-  C.expression = rhs;
-  ramsey_constraints.push_back(C);
+  reset_data_tree();
 }
 
 void
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index a664d1c5..0c4fe58c 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -186,8 +186,6 @@ private:
   MomentCalibration::constraints_t moment_calibration_constraints;
   //! Temporary storage for irf_calibration
   IrfCalibration::constraints_t irf_calibration_constraints;
-  //! Temporary storage for ramsey_constraints
-  RamseyConstraintsStatement::constraints_t ramsey_constraints;
   //! Temporary storage for svar_identification blocks
   SvarIdentificationStatement::svar_identification_restrictions_t svar_ident_restrictions;
   //! Temporary storage for mapping the equation number to the restrictions within an
@@ -667,18 +665,9 @@ public:
   void end_planner_objective(expr_t expr);
   //! Ramsey model statement
   void ramsey_model();
-  //! Ramsey constraints statement
-  void add_ramsey_constraints_statement();
-  //! Ramsey less constraint
-  void ramsey_constraint_add_less(const string& name, const expr_t rhs);
-  //! Ramsey greater constraint
-  void ramsey_constraint_add_greater(const string& name, const expr_t rhs);
-  //! Ramsey less or equal constraint
-  void ramsey_constraint_add_less_equal(const string& name, const expr_t rhs);
-  //! Ramsey greater or equal constraint
-  void ramsey_constraint_add_greater_equal(const string& name, const expr_t rhs);
-  //! Ramsey constraint helper function
-  void add_ramsey_constraint(const string& name, BinaryOpcode op_code, const expr_t rhs);
+  // Ramsey constraints statement
+  void begin_ramsey_constraints();
+  void end_ramsey_constraints(const vector<expr_t>& constraints);
   //! Ramsey policy statement
   void ramsey_policy(vector<string> symbol_list);
   //! Evaluate Planner Objective
diff --git a/src/Statement.hh b/src/Statement.hh
index df6cd1f0..e67525cb 100644
--- a/src/Statement.hh
+++ b/src/Statement.hh
@@ -169,8 +169,6 @@ struct ModFileStructure
   bool endval_learnt_in_present {false};
   // Whether an occbin_constraints block appears
   bool occbin_constraints_present {false};
-  // Whether a ramsey_constraints block appears
-  bool ramsey_constraints_present {false};
 };
 
 class Statement
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index a79a1aa9..a1d68819 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -87,17 +87,18 @@ StaticModel::StaticModel(const DynamicModel& m) :
         if (dynamic_equations.contains(i))
           {
             auto [static_only_equations, static_only_equations_lineno,
-                  static_only_equations_equation_tags]
+                  static_only_complementarity_conditions, static_only_equations_equation_tags]
                 = m.getStaticOnlyEquationsInfo();
 
             addEquation(static_only_equations[static_only_index]->toStatic(*this),
                         static_only_equations_lineno[static_only_index],
+                        static_only_complementarity_conditions[static_only_index],
                         static_only_equations_equation_tags.getTagsByEqn(static_only_index));
             static_only_index++;
           }
         else
           addEquation(m.equations[i]->toStatic(*this), m.equations_lineno[i],
-                      m.equation_tags.getTagsByEqn(i));
+                      m.complementarity_conditions[i], m.equation_tags.getTagsByEqn(i));
       }
     catch (DataTree::DivisionByZeroException)
       {
@@ -250,6 +251,8 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder,
            << endl;
       exit(EXIT_FAILURE);
     }
+
+  computeMCPEquationsReordering();
 }
 
 void
@@ -305,6 +308,8 @@ StaticModel::writeStaticFile(const string& basename, bool use_dll, const string&
 
   writeSetAuxiliaryVariablesFile<false>(basename, julia);
 
+  writeComplementarityConditionsFile<false>(basename);
+
   // Support for model debugging
   if (!julia)
     writeDebugModelMFiles<false>(basename);
@@ -331,6 +336,11 @@ StaticModel::writeDriverOutput(ostream& output) const
     writeBlockDriverOutput(output);
 
   writeDriverSparseIndicesHelper<false, false>(output);
+
+  output << "M_.static_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
 }
 
 void
-- 
GitLab