#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2023 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
#################################################################################
"""
IDAES Homotopy meta-solver routine.
"""
__author__ = "Andrew Lee"
import logging
from pyomo.environ import Block, SolverFactory, TerminationCondition
from pyomo.core.base.var import _VarData
from pyomo.contrib.parmest.utils.ipopt_solver_wrapper import ipopt_solve_with_stats
from idaes.core.util.model_serializer import to_json, from_json
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.util.exceptions import ConfigurationError
import idaes.logger as idaeslog
_log = idaeslog.getLogger(__name__)
[docs]def homotopy(
model,
variables,
targets,
max_solver_iterations=50,
max_solver_time=10,
step_init=0.1,
step_cut=0.5,
iter_target=4,
step_accel=0.5,
max_step=1,
min_step=0.05,
max_eval=200,
):
"""
Homotopy meta-solver routine using Ipopt as the non-linear solver. This
routine takes a model along with a list of fixed variables in that model
and a list of target values for those variables. The routine then tries to
iteratively move the values of the fixed variables to their target values
using an adaptive step size.
Args:
model : model to be solved
variables : list of Pyomo Var objects to be varied using homotopy.
Variables must be fixed.
targets : list of target values for each variable
max_solver_iterations : maximum number of solver iterations per
homotopy step (default=50)
max_solver_time : maximum cpu time for the solver per homotopy step
(default=10)
step_init : initial homotopy step size (default=0.1)
step_cut : factor by which to reduce step size on failed step
(default=0.5)
step_accel : acceleration factor for adjusting step size on successful
step (default=0.5)
iter_target : target number of solver iterations per homotopy step
(default=4)
max_step : maximum homotopy step size (default=1)
min_step : minimum homotopy step size (default=0.05)
max_eval : maximum number of homotopy evaluations (both successful and
unsuccessful) (default=200)
Returns:
Termination Condition : A Pyomo TerminationCondition Enum indicating
how the meta-solver terminated (see documentation)
Solver Progress : a fraction indication how far the solver progressed
from the initial values to the target values
Number of Iterations : number of homotopy evaluations before solver
terminated
"""
eps = 1e-3 # Tolerance for homotopy step convergence to 1
# Get model logger
_log = logging.getLogger(__name__)
# Validate model is an instance of Block
if not isinstance(model, Block):
raise TypeError(
"Model provided was not a valid Pyomo model object "
"(instance of Block). Please provide a valid model."
)
if degrees_of_freedom(model) != 0:
raise ConfigurationError(
"Degrees of freedom in model are not equal to zero. Homotopy "
"should not be used on problems which are not well-defined."
)
# Validate variables and targets
if len(variables) != len(targets):
raise ConfigurationError("Number of variables and targets do not match.")
for i, v in enumerate(variables):
t = targets[i]
if not isinstance(v, _VarData):
raise TypeError(
"Variable provided ({}) was not a valid Pyomo Var "
"component.".format(v)
)
# Check that v is part of model
parent = v.parent_block()
while parent != model:
if parent is None:
raise ConfigurationError("Variable {} is not part of model".format(v))
parent = parent.parent_block()
# Check that v is fixed
if not v.fixed:
raise ConfigurationError(
"Homotopy metasolver provided with unfixed variable {}."
"All variables must be fixed.".format(v.name)
)
# Check bounds on v (they don't really matter, but check for sanity)
if v.ub is not None:
if v.value > v.ub:
raise ConfigurationError(
"Current value for variable {} is greater than the "
"upper bound for that variable. Please correct this "
"before continuing.".format(v.name)
)
if t > v.ub:
raise ConfigurationError(
"Target value for variable {} is greater than the "
"upper bound for that variable. Please correct this "
"before continuing.".format(v.name)
)
if v.lb is not None:
if v.value < v.lb:
raise ConfigurationError(
"Current value for variable {} is less than the "
"lower bound for that variable. Please correct this "
"before continuing.".format(v.name)
)
if t < v.lb:
raise ConfigurationError(
"Target value for variable {} is less than the "
"lower bound for that variable. Please correct this "
"before continuing.".format(v.name)
)
# TODO : Should we be more restrictive on these values to avoid users
# TODO : picking numbers that are less likely to solve (but still valid)?
# Validate homotopy parameter selections
if not 0.05 <= step_init <= 0.8:
raise ConfigurationError(
"Invalid value for step_init ({}). Must lie "
"between 0.05 and 0.8.".format(step_init)
)
if not 0.1 <= step_cut <= 0.9:
raise ConfigurationError(
"Invalid value for step_cut ({}). Must lie "
"between 0.1 and 0.9.".format(step_cut)
)
if step_accel < 0:
raise ConfigurationError(
"Invalid value for step_accel ({}). Must be "
"greater than or equal to 0.".format(step_accel)
)
if iter_target < 1:
raise ConfigurationError(
"Invalid value for iter_target ({}). Must be "
"greater than or equal to 1.".format(iter_target)
)
if not isinstance(iter_target, int):
raise ConfigurationError(
"Invalid value for iter_target ({}). Must be "
"an an integer.".format(iter_target)
)
if not 0.05 <= max_step <= 1:
raise ConfigurationError(
"Invalid value for max_step ({}). Must lie "
"between 0.05 and 1.".format(max_step)
)
if not 0.01 <= min_step <= 0.1:
raise ConfigurationError(
"Invalid value for min_step ({}). Must lie "
"between 0.01 and 0.1.".format(min_step)
)
if min_step > max_step:
raise ConfigurationError(
"Invalid arguments: step_min must be less or equal to step_max."
)
if not min_step <= step_init <= max_step:
raise ConfigurationError(
"Invalid arguments: step_init must lie between min_step and max_step."
)
if max_eval < 1:
raise ConfigurationError(
"Invalid value for max_eval ({}). Must be "
"greater than or equal to 1.".format(step_accel)
)
if not isinstance(max_eval, int):
raise ConfigurationError(
"Invalid value for max_eval ({}). Must be "
"an an integer.".format(iter_target)
)
# Create solver object
solver_obj = SolverFactory("ipopt")
# Perform initial solve of model to confirm feasible initial solution
(
results, # pylint: disable=unused-variable
solved,
sol_iter,
sol_time, # pylint: disable=unused-variable
sol_reg,
) = ipopt_solve_with_stats(
model, solver_obj, max_solver_iterations, max_solver_time
)
if not solved:
_log.exception("Homotopy Failed - initial solution infeasible.")
return TerminationCondition.infeasible, 0, 0
elif sol_reg != "-":
_log.warning("Homotopy - initial solution converged with regularization.")
return TerminationCondition.other, 0, 0
else:
_log.info("Homotopy - initial point converged")
# Set up homotopy variables
# Get initial values and deltas for all variables
v_init = []
for i, v in enumerate(variables):
v_init.append(v.value)
n_0 = 0.0 # Homotopy progress variable
s = step_init # Set step size to step_init
iter_count = 0 # Counter for homotopy iterations
# Save model state to dict
# TODO : for very large models, it may be necessary to dump this to a file
current_state = to_json(model, return_dict=True)
while n_0 < 1.0:
iter_count += 1 # Increase iter_count regardless of success or failure
# Calculate next n value given current step size
if n_0 + s >= 1.0 - eps:
n_1 = 1.0
else:
n_1 = n_0 + s
_log.info(f"Homotopy Iteration {iter_count}. Next Step: {n_1} (Current: {n_0})")
# Update values for all variables using n_1
for i, v in enumerate(variables):
v.fix(targets[i] * n_1 + v_init[i] * (1 - n_1))
# Solve model at new state
results, solved, sol_iter, sol_time, sol_reg = ipopt_solve_with_stats(
model, solver_obj, max_solver_iterations, max_solver_time
)
# Check solver output for convergence
if solved:
# Step succeeded - accept current state
current_state = to_json(model, return_dict=True)
# Update n_0 to accept current step
n_0 = n_1
# Check solver iterations and calculate next step size
s_proposed = s * (1 + step_accel * (iter_target / sol_iter - 1))
if s_proposed > max_step:
s = max_step
elif s_proposed < min_step:
s = min_step
else:
s = s_proposed
else:
# Step failed - reload old state
from_json(model, current_state)
# Try to cut back step size
if s > min_step:
# Step size can be cut
s = max(min_step, s * step_cut)
else:
# Step is already at minimum size, terminate homotopy
_log.exception(
f"Homotopy failed - could not converge at minimum step "
f"size. Current progress is {n_0}"
)
return TerminationCondition.minStepLength, n_0, iter_count
if iter_count >= max_eval: # Use greater than or equal to to be safe
_log.exception(
f"Homotopy failed - maximum homotopy iterations "
f"exceeded. Current progress is {n_0}"
)
return TerminationCondition.maxEvaluations, n_0, iter_count
if sol_reg == "-":
_log.info(
f"Homotopy successful - converged at target values in {iter_count} iterations."
)
return TerminationCondition.optimal, n_0, iter_count
else:
_log.exception(
f"Homotopy failed - converged at target values with "
f"regularization in {iter_count} iterations."
)
return TerminationCondition.other, n_0, iter_count
