Solving Dynamic Optimization Problems with Schur-Complement Decomposition

In order to solve a dynamic optimization problem with schur-complement decomposition, you must create a class which inherits from MPIDynamicSchurComplementInteriorPointInterface. This class must implement the method build_model_for_time_block():

import parapint

class Problem(parapint.interfaces.MPIDynamicSchurComplementInteriorPointInterface):
    def __init__(self, your_arguments):
        # do anything you need to here
        super(Problem, self).__init__(start_t, end_t, num_time_blocks, mpi_comm)

    def build_model_for_time_block(self, ndx, start_t, end_t, add_init_conditions):
        # build the dynamic optimization problem with Pyomo over the time horizon
        # [start_t, end_t] and return the model along with two lists. The first
        # list should be a list of pyomo variables corresponding to the states at
        # start_t. The second list should be a list of pyomo variables
        # corresponding to the states at end_t. Continuity will be enforced
        # between the states at end_t for one time block
        # and the states at start_t for the next time block. It is very important for
        # the ordering of the state variables to be the same for every time block.

        return model, start_states, end_states

problem = Problem(some_arguments)

The build_model_for_time_block() method will be called once for every time block. It will be called within the call to __init__() of the super class (MPIDynamicSchurComplementInteriorPointInterface). Therefore, if you override the __init__ method, it is very important to still call the __init__ method of the base class as shown above. There is an example class in schur_complement.py in the examples directory within Parapint.

In addition to the implementation of the class described above, you must create an instance of MPISchurComplementLinearSolver. This linear solver requires a sub-solver for every time block:

cntl_options = {1: 1e-6}  # the pivot tolerance
sub_solvers = {ndx: parapint.linalg.InteriorPointMA27Interface(cntl_options=cntl_options) for ndx in range(num_time_blocks)}
schur_complement_solver = parapint.linalg.InteriorPointMA27Interface(cntl_options=cntl_options)
linear_solver = parapint.linalg.MPISchurComplementLinearSolver(subproblem_solvers=sub_solvers,
                                                              schur_complement_solver=schur_complement_solver)

The linear solver and interface instances can then be used with the interior point algorithm:

options = parapint.algorithms.IPOptions()
options.linalg.solver = linear_solver
status = parapint.algorithms.ip_solve(interface, options)
assert status == parapint.interior_point.InteriorPointStatus.optimal
problem.load_primals_into_pyomo_model()
for ndx in problem.local_block_indices:
    model = problem.pyomo_model(ndx)
    model.pprint()