diff --git a/pyomo/contrib/gdpopt/tests/four_stage_dynamic_model.py b/pyomo/contrib/gdpopt/tests/four_stage_dynamic_model.py new file mode 100644 index 00000000000..2d4dbe9c02c --- /dev/null +++ b/pyomo/contrib/gdpopt/tests/four_stage_dynamic_model.py @@ -0,0 +1,387 @@ +# from pyomo.environ import * +from pyomo.core import ( + Var, + Constraint, + Objective, + Set, + minimize, + exp, + ConcreteModel, + LogicalConstraint, + exactly, + lnot, + lor, + BooleanVar, + land, +) +from pyomo.dae import Integral, DerivativeVar, ContinuousSet +from pyomo.gdp import Disjunct, Disjunction + + +def build_model(mode_transfer=False): + model = ConcreteModel() + + # Set + model.stage = Set(initialize=[1, 2, 3, 4]) + model.mode = Set(initialize=[1, 2, 3]) + + model.t1 = ContinuousSet(bounds=(0, 1)) + model.t2 = ContinuousSet(bounds=(1, 2)) + model.t3 = ContinuousSet(bounds=(2, 3)) + model.t4 = ContinuousSet(bounds=(3, 4)) + + # Variables + model.x1 = Var(model.t1, bounds=(0, 10)) + model.x2 = Var(model.t2, bounds=(0, 10)) + model.x3 = Var(model.t3, bounds=(0, 10)) + model.x4 = Var(model.t4, bounds=(0, 10)) + model.u1 = Var(bounds=(-4, 4)) + model.u2 = Var(bounds=(-4, 4)) + model.u3 = Var(bounds=(-4, 4)) + model.u4 = Var(bounds=(-4, 4)) + + # Dynamic model + model.dxdt1 = DerivativeVar(model.x1, wrt=model.t1) + model.dxdt2 = DerivativeVar(model.x2, wrt=model.t2) + model.dxdt3 = DerivativeVar(model.x3, wrt=model.t3) + model.dxdt4 = DerivativeVar(model.x4, wrt=model.t4) + + # logic constraint + model.stage_mode = Disjunct(model.stage * model.mode) + model.d = Disjunction(model.stage) + + # Stage 1 + + def stage1_mode1_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt1[t] == -model.x1[t] * exp(model.x1[t] - 1) + model.u1 + + model.stage_mode[1, 1].mode1_dynamic_constraint = Constraint( + model.t1, rule=stage1_mode1_dynamic + ) + + def stage1_mode2_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt1[t] == (0.5 * model.x1[t] ** 3 + model.u1) / 20 + + model.stage_mode[1, 2].mode2_dynamic_constraint = Constraint( + model.t1, rule=stage1_mode2_dynamic + ) + + def stage1_mode3_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt1[t] == (model.x1[t] ** 2 + model.u1) / (t + 20) + + model.stage_mode[1, 3].mode3_dynamic_constraint = Constraint( + model.t1, rule=stage1_mode3_dynamic + ) + + # Stage 2 + + def stage2_mode1_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt2[t] == -model.x2[t] * exp(model.x2[t] - 1) + model.u2 + + model.stage_mode[2, 1].mode1_dynamic_constraint = Constraint( + model.t2, rule=stage2_mode1_dynamic + ) + + def stage2_mode2_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt2[t] == (0.5 * model.x2[t] ** 3 + model.u2) / 20 + + model.stage_mode[2, 2].mode2_dynamic_constraint = Constraint( + model.t2, rule=stage2_mode2_dynamic + ) + + def stage2_mode3_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt2[t] == (model.x2[t] ** 2 + model.u2) / (t + 20) + + model.stage_mode[2, 3].mode3_dynamic_constraint = Constraint( + model.t2, rule=stage2_mode3_dynamic + ) + + # Stage 3 + + def stage3_mode1_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt3[t] == -model.x3[t] * exp(model.x3[t] - 1) + model.u3 + + model.stage_mode[3, 1].mode1_dynamic_constraint = Constraint( + model.t3, rule=stage3_mode1_dynamic + ) + + def stage3_mode2_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt3[t] == (0.5 * model.x3[t] ** 3 + model.u3) / 20 + + model.stage_mode[3, 2].mode2_dynamic_constraint = Constraint( + model.t3, rule=stage3_mode2_dynamic + ) + + def stage3_mode3_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt3[t] == (model.x3[t] ** 2 + model.u3) / (t + 20) + + model.stage_mode[3, 3].mode3_dynamic_constraint = Constraint( + model.t3, rule=stage3_mode3_dynamic + ) + + # Stage 4 + + def stage4_mode1_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt4[t] == -model.x4[t] * exp(model.x4[t] - 1) + model.u4 + + model.stage_mode[4, 1].mode1_dynamic_constraint = Constraint( + model.t4, rule=stage4_mode1_dynamic + ) + + def stage4_mode2_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt4[t] == (0.5 * model.x4[t] ** 3 + model.u4) / 20 + + model.stage_mode[4, 2].mode2_dynamic_constraint = Constraint( + model.t4, rule=stage4_mode2_dynamic + ) + + def stage4_mode3_dynamic(disjunct, t): + model = disjunct.model() + return model.dxdt4[t] == (model.x4[t] ** 2 + model.u4) / (t + 20) + + model.stage_mode[4, 3].mode3_dynamic_constraint = Constraint( + model.t4, rule=stage4_mode3_dynamic + ) + + model.d[1] = [ + model.stage_mode[1, 1], + model.stage_mode[1, 2], + model.stage_mode[1, 3], + ] + model.d[2] = [ + model.stage_mode[2, 1], + model.stage_mode[2, 2], + model.stage_mode[2, 3], + ] + model.d[3] = [ + model.stage_mode[3, 1], + model.stage_mode[3, 2], + model.stage_mode[3, 3], + ] + model.d[4] = [ + model.stage_mode[4, 1], + model.stage_mode[4, 2], + model.stage_mode[4, 3], + ] + + if mode_transfer: + model.lc1 = LogicalConstraint( + expr=exactly( + 1, + model.stage_mode[1, 1].indicator_var, + model.stage_mode[1, 2].indicator_var, + model.stage_mode[1, 3].indicator_var, + ) + ) + model.lc2 = LogicalConstraint( + expr=exactly( + 1, + model.stage_mode[2, 1].indicator_var, + model.stage_mode[2, 2].indicator_var, + model.stage_mode[2, 3].indicator_var, + ) + ) + model.lc3 = LogicalConstraint( + expr=exactly( + 1, + model.stage_mode[3, 1].indicator_var, + model.stage_mode[3, 2].indicator_var, + model.stage_mode[3, 3].indicator_var, + ) + ) + model.lc4 = LogicalConstraint( + expr=exactly( + 1, + model.stage_mode[4, 1].indicator_var, + model.stage_mode[4, 2].indicator_var, + model.stage_mode[4, 3].indicator_var, + ) + ) + model.transfer_stage1 = Set(initialize=[2, 3, 4]) + model.transfer_stage2 = Set(initialize=[2, 3, 4, 5]) + model.mode_stransfer_set = Set(initialize=[1, 2]) + model.mode_transfer = BooleanVar( + model.transfer_stage2, model.mode_stransfer_set + ) + model.mode_transfer_lc1 = LogicalConstraint( + expr=exactly( + 1, + model.mode_transfer[2, 1], + model.mode_transfer[3, 1], + model.mode_transfer[4, 1], + model.mode_transfer[5, 1], + ) + ) + model.mode_transfer_lc2 = LogicalConstraint( + expr=exactly( + 1, + model.mode_transfer[2, 2], + model.mode_transfer[3, 2], + model.mode_transfer[4, 2], + model.mode_transfer[5, 2], + ) + ) + + def _mode_transfer_rule1(model, stage): + return model.mode_transfer[stage, 1].equivalent_to( + land( + model.stage_mode[stage - 1, 1].indicator_var, + model.stage_mode[stage, 2].indicator_var, + ) + ) + + model.mode_transfer2mode_choice_lc1 = LogicalConstraint( + model.transfer_stage1, rule=_mode_transfer_rule1 + ) + + def _mode_transfer_rule2(model, stage): + return model.mode_transfer[stage, 2].equivalent_to( + land( + model.stage_mode[stage - 1, 2].indicator_var, + model.stage_mode[stage, 3].indicator_var, + ) + ) + + model.mode_transfer2mode_choice_lc2 = LogicalConstraint( + model.transfer_stage1, rule=_mode_transfer_rule2 + ) + + def _mode_transfer_rule3(model, stage): + return model.mode_transfer[stage, 2].implies( + lor( + model.mode_transfer[stage1, 1] + for stage1 in model.transfer_stage1 + if stage1 < stage + ) + ) + + model.mode_transfer2mode_choice_lc3 = LogicalConstraint( + model.transfer_stage1, rule=_mode_transfer_rule3 + ) + + def _mode_transfer_rule4(model): + return model.mode_transfer[5, 1].implies( + lnot( + lor( + model.stage_mode[stage1, 2].indicator_var + for stage1 in model.stage + ) + ) + ) + + model.mode_transfer2mode_choice_lc4 = LogicalConstraint( + rule=_mode_transfer_rule4 + ) + + def _mode_transfer_rule5(model): + return model.mode_transfer[5, 2].implies( + lnot( + lor( + model.stage_mode[stage1, 3].indicator_var + for stage1 in model.stage + ) + ) + ) + + model.mode_transfer2mode_choice_lc5 = LogicalConstraint( + rule=_mode_transfer_rule5 + ) + + # Sequence constraint + def _sequence_rule1(model, stage): + if stage == 1: + return Constraint.Skip + else: + return model.stage_mode[stage, 2].indicator_var.implies( + lor( + model.stage_mode[stage2, 1].indicator_var + for stage2 in model.stage + if stage2 < stage + ) + ) + + model.seq1 = LogicalConstraint(model.stage, rule=_sequence_rule1) + model.stage_mode[1, 2].indicator_var.fix(False) + + def _sequence_rule2(model, stage): + if stage == 4: + return Constraint.Skip + else: + return model.stage_mode[stage, 2].indicator_var.implies( + lnot( + lor( + model.stage_mode[stage2, 1].indicator_var + for stage2 in model.stage + if stage2 > stage + ) + ) + ) + + model.seq2 = LogicalConstraint(model.stage, rule=_sequence_rule2) + + def _sequence_rule3(model, stage): + if stage <= 1: + return Constraint.Skip + else: + return model.stage_mode[stage, 3].indicator_var.implies( + lor( + model.stage_mode[stage2, 2].indicator_var + for stage2 in model.stage + if stage2 < stage + ) + ) + + model.seq3 = LogicalConstraint(model.stage, rule=_sequence_rule3) + model.stage_mode[1, 3].indicator_var.fix(False) + model.stage_mode[2, 3].indicator_var.fix(False) + + def _sequence_rule4(model, stage): + if stage == 4: + return Constraint.Skip + else: + return model.stage_mode[stage, 3].indicator_var.implies( + lnot( + lor( + model.stage_mode[stage2, 2].indicator_var + for stage2 in model.stage + if stage2 > stage + ) + ) + ) + + model.seq4 = LogicalConstraint(model.stage, rule=_sequence_rule4) + + model.c1 = Constraint(expr=model.x1[0] == 1) + model.c2 = Constraint(expr=model.x1[1] == model.x2[1]) + model.c3 = Constraint(expr=model.x2[2] == model.x3[2]) + model.c4 = Constraint(expr=model.x3[3] == model.x4[3]) + + # Objective function + model.intx1 = Integral( + model.t1, wrt=model.t1, rule=lambda model, t: model.x1[t] ** 2 + ) + model.intx2 = Integral( + model.t2, wrt=model.t2, rule=lambda model, t: model.x2[t] ** 2 + ) + model.intx3 = Integral( + model.t3, wrt=model.t3, rule=lambda model, t: model.x3[t] ** 2 + ) + model.intx4 = Integral( + model.t4, wrt=model.t4, rule=lambda model, t: model.x4[t] ** 2 + ) + + model.obj = Objective( + expr=-(model.intx1 + model.intx2 + model.intx3 + model.intx4), sense=minimize + ) + return model diff --git a/pyomo/contrib/gdpopt/tests/test_ldsda.py b/pyomo/contrib/gdpopt/tests/test_ldsda.py new file mode 100644 index 00000000000..9bcb2fd10cb --- /dev/null +++ b/pyomo/contrib/gdpopt/tests/test_ldsda.py @@ -0,0 +1,47 @@ +from pyomo.environ import SolverFactory, value, Var, Constraint, TransformationFactory +from pyomo.gdp import Disjunct +import pyomo.common.unittest as unittest +from four_stage_dynamic_model import build_model + + +class TestGDPoptLDSDA(unittest.TestCase): + """Real unit tests for GDPopt""" + + @unittest.skipUnless(SolverFactory('gams').available(), "gams solver not available") + def test_solve_four_stage_dynamic_model(self): + + model = build_model(mode_transfer=True) + + # Discretize the model using dae.collocation + discretizer = TransformationFactory('dae.collocation') + discretizer.apply_to(model, nfe=10, ncp=3, scheme='LAGRANGE-RADAU') + # We need to reconstruct the constraints in disjuncts after discretization. + # This is a bug in Pyomo.dae. https://github.com/Pyomo/pyomo/issues/3101 + for disjunct in model.component_data_objects(ctype=Disjunct): + for constraint in disjunct.component_objects(ctype=Constraint): + constraint._constructed = False + constraint.construct() + + for dxdt in model.component_data_objects(ctype=Var, descend_into=True): + if 'dxdt' in dxdt.name: + dxdt.setlb(-300) + dxdt.setub(300) + + for direction_norm in ['L2', 'Linf']: + result = SolverFactory('gdpopt.ldsda').solve( + model, + direction_norm=direction_norm, + minlp_solver='gams', + minlp_solver_args=dict(solver='knitro'), + starting_point=[1, 2], + logical_constraint_list=[ + model.mode_transfer_lc1.name, + model.mode_transfer_lc2.name, + ], + time_limit=100, + ) + self.assertAlmostEqual(value(model.obj), -23.305325, places=4) + + +if __name__ == '__main__': + unittest.main()