diff --git a/svi/auto_thermal_reformer/compute_surrogate_error.py b/svi/auto_thermal_reformer/compute_surrogate_error.py new file mode 100644 index 0000000..630258f --- /dev/null +++ b/svi/auto_thermal_reformer/compute_surrogate_error.py @@ -0,0 +1,127 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import pyomo.environ as pyo +from pyomo.common.timing import TicTocTimer +from pyomo.util.subsystems import TemporarySubsystemManager +from pyomo.contrib.incidence_analysis import solve_strongly_connected_components +from idaes.core.util.exceptions import InitializationError +from svi.auto_thermal_reformer.reactor_model import create_instance + + +def compute_surrogate_error(model): + model_surrogate = model.fs.reformer_surrogate + + model_surrogate_outputs = { + key: var.value for key, var in model_surrogate.output_vars_as_dict().items() + } + + inputs = [ + model.fs.reformer_bypass.reformer_outlet.flow_mol[0], + model.fs.reformer_bypass.reformer_outlet.temperature[0], + model.fs.reformer_mix.steam_inlet.flow_mol[0], + model.fs.reformer.conversion, + ] + + scc_solver = pyo.SolverFactory("ipopt") + + timer = TicTocTimer() + timer.tic() + + with TemporarySubsystemManager(to_fix=inputs): + surrogate_copy = model_surrogate.clone() + timer.toc("clone-surrogate") + solve_strongly_connected_components(surrogate_copy, solver=scc_solver) + timer.toc("solve-scc-surrogate") + surrogate_res = scc_solver.solve(surrogate_copy) + if not pyo.check_optimal_termination(surrogate_res): + print("WARNING: Surrogate model failed to simulate") + raise ValueError("Surrogate reactor model failed to simulate") + + # If we have converged, these newly computed surrogate outputs should be + # the same as our "model surrogate outputs" above. If we did not converge, + # they may be different. + surrogate_outputs = { + key: var.value for key, var in surrogate_copy.output_vars_as_dict().items() + } + + try: + fullspace_model = create_instance( + model.fs.reformer.conversion.value, + model.fs.reformer_mix.steam_inlet.flow_mol[0].value, + # Note that "reformer_outlet" here means "the outlet of the bypass splitter + # that goes to the reformer". The other outlet is called "bypass_outlet". + # The inlet is simply called "inlet". + model.fs.reformer_bypass.reformer_outlet.flow_mol[0].value, + model.fs.reformer_bypass.reformer_outlet.temperature[0].value, + initialize=True, + ) + timer.toc("create-instance-fullspace") + except InitializationError: + print("WARNING: Full-space model failed to initialize. Trying to continue.") + fullspace_model = create_instance( + model.fs.reformer.conversion.value, + model.fs.reformer_mix.steam_inlet.flow_mol[0].value, + model.fs.reformer_bypass.reformer_outlet.flow_mol[0].value, + model.fs.reformer_bypass.reformer_outlet.temperature[0].value, + initialize=False, + ) + + solve_strongly_connected_components( + fullspace_model, + solver=scc_solver, + use_calc_var=False, + ) + timer.toc("solve-scc-fullspace") + fullspace_res = scc_solver.solve(fullspace_model, tee=True) + timer.toc("solve-fullspace") + + if not pyo.check_optimal_termination(fullspace_res): + # In parameter sweep scripts, this ValueError is caught and we + # write an empty row in the surrogate-error file. + # But we don't want to count a failure for the surrogate if this + # fails... this can be handled by the caller. + print("WARNING: Full-space model failed to simulate") + raise ValueError("Full-space reactor model failed to simulate") + + fullspace_output = { + "Fout": fullspace_model.fs.reformer.outlet.flow_mol[0].value, + "Tout": fullspace_model.fs.reformer.outlet.temperature[0].value, + "HeatDuty": fullspace_model.fs.reformer.heat_duty[0].value, + } + for key in surrogate_outputs: + # The keys we have not added so far are component names, which may + # be used as indices to mole_frac_comp + if key not in fullspace_output: + fullspace_output[key] = fullspace_model.fs.reformer.outlet.mole_frac_comp[0, key].value + + relative_errors = { + key: ( + abs(fullspace_output[key] - surrogate_outputs[key]) + / max(1, fullspace_output[key], surrogate_outputs[key]) + ) + for key in surrogate_outputs + } + #max_relative_error = max(relative_errors.values()) + #ave_relative_error = sum(relative_errors.values()) / len(relative_errors) + timer.toc("done") + + return relative_errors diff --git a/svi/auto_thermal_reformer/condition_num_plots.py b/svi/auto_thermal_reformer/condition_num_plots.py new file mode 100644 index 0000000..59bd3de --- /dev/null +++ b/svi/auto_thermal_reformer/condition_num_plots.py @@ -0,0 +1,194 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import pyomo.environ as pyo +from svi.cyipopt import ConditioningCallback +from pyomo.contrib.incidence_analysis import solve_strongly_connected_components +from svi.auto_thermal_reformer.fullspace_flowsheet import ( + make_optimization_model, + make_simulation_model, +) +from svi.auto_thermal_reformer.implicit_flowsheet import make_implicit +from svi.external import add_external_function_libraries_to_environment +from pyomo.contrib.pynumero.interfaces.external_pyomo_model import ExternalPyomoModel +from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock +from pyomo.contrib.pynumero.algorithms.solvers.implicit_functions import ( + CyIpoptSolverWrapper +) + +from svi.auto_thermal_reformer.alamo_flowsheet import ( + create_instance as create_alamo_instance, + initialize_alamo_atr_flowsheet, + DEFAULT_SURROGATE_FNAME as DEFAULT_ALAMO_SURROGATE_FNAME +) + +from svi.auto_thermal_reformer.nn_flowsheet import ( + create_instance as create_nn_instance, + initialize_nn_atr_flowsheet, + DEFAULT_SURROGATE_FNAME as DEFAULT_NN_SURROGATE_FNAME, +) +import svi.auto_thermal_reformer.config as config +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import argparse + +# SUBSETS represent a tuple containing (X,P,iters), where iters = number of iterations +# for the {full-space, implicit function, ALAMO, NN} formulations to successfully converge. +# For the unsuccessful instances, implicit function displays eval error after iter 12. The +# full space and NN formulations do not find a solution after 300 iters. + +data_dir = config.get_data_dir() + +SUCCESS_SUBSETS_FULLSPACE = [ + (0.94, 1947379.0, 89), + (0.90, 1727379.0, 123), + (0.92, 1587379.0, 58), + (0.93, 1657379.0, 52) +] + +SUCCESS_SUBSETS_IMPLICIT = [ + (0.94, 1947379.0, 45), + (0.90, 1727379.0, 35), + (0.92, 1587379.0, 41), + (0.93, 1657379.0, 37) +] + +SUCCESS_SUBSETS_ALAMO = [ + (0.94, 1947379.0, 29), + (0.90, 1727379.0, 26), + (0.92, 1587379.0, 32), + (0.93, 1657379.0, 27) +] + +SUCCESS_SUBSETS_NN = [ + (0.94, 1947379.0, 53), + (0.90, 1727379.0, 48), + (0.92, 1587379.0, 47), + (0.93, 1657379.0, 44) +] + +NOSUCCESS_SUBSETS_FULLSPACE = [ + (0.96, 1447379.0, 150), + (0.97, 1447379.0, 150), + (0.96, 1587379.0, 150), + (0.97, 1587379.0, 150) +] + +NOSUCCESS_SUBSETS_IMPLICIT = [ + (0.96, 1447379.0, 12), + (0.97, 1447379.0, 12), + (0.96, 1587379.0, 12), + (0.97, 1587379.0, 12) +] + +NOSUCCESS_SUBSETS_NN = [ + (0.96, 1447379.0, 150), + (0.97, 1447379.0, 150), + (0.96, 1587379.0, 150), + (0.97, 1587379.0, 150) +] + +def solver(m, iters = 300): + solver = config.get_optimization_solver(iters = iters) + callback = ConditioningCallback() + solver.config.intermediate_callback = callback + solver.solve(m, tee=True) + condition_numbers_list = callback.condition_numbers + return condition_numbers_list[1:] + +def full(X=0.94, P=1947379.0, iters=300): + m = make_optimization_model(X, P) + print(f"Solving sample with X={X}, P={P}") + condition_numbers_list = solver(m, iters = iters) + return condition_numbers_list + +def implicit(X=0.94, P=1947379.0, iters=300): + m = make_optimization_model(X, P) + add_external_function_libraries_to_environment(m) + m_implicit = make_implicit(m) + print(f"Solving sample with X={X}, P={P}") + condition_numbers_list = solver(m_implicit, iters = iters) + return condition_numbers_list + +def alamo(X=0.94, P=1947379.0, iters=300): + m = create_alamo_instance(X, P, surrogate_fname=os.path.join(data_dir, DEFAULT_ALAMO_SURROGATE_FNAME)) + initialize_alamo_atr_flowsheet(m) + m.fs.reformer_bypass.inlet.temperature.unfix() + m.fs.reformer_bypass.inlet.flow_mol.unfix() + print(f"Solving sample with X={X}, P={P}") + condition_numbers_list = solver(m, iters = iters) + return condition_numbers_list + +def nn(X=0.94, P=1947379.0, iters=300): + m = create_nn_instance(X, P, surrogate_fname=os.path.join(data_dir, DEFAULT_NN_SURROGATE_FNAME)) + initialize_nn_atr_flowsheet(m) + m.fs.reformer_bypass.inlet.temperature.unfix() + m.fs.reformer_bypass.inlet.flow_mol.unfix() + print(f"Solving sample with X={X}, P={P}") + condition_numbers_list = solver(m, iters = iters) + return condition_numbers_list + +def plot_subsets(SUBSETS_FULLSPACE, SUBSETS_IMPLICIT, SUBSETS_NN, SUBSETS_ALAMO=None, unsuccessful=False): + plt.figure(figsize=(12, 8)) + + for i, (full_subset, implicit_subset, nn_subset) in enumerate(zip(SUBSETS_FULLSPACE, SUBSETS_IMPLICIT, SUBSETS_NN), 1): + plt.subplot(2, 2, i) + full_list = full(X=full_subset[0], P=full_subset[1], iters=full_subset[2]) + implicit_list = implicit(X=implicit_subset[0], P=implicit_subset[1], iters=implicit_subset[2]) + nn_list = nn(X=nn_subset[0], P=nn_subset[1], iters=nn_subset[2]) + plt.plot(full_list, label='Full-space') + plt.plot(implicit_list, label='Implicit function') + plt.plot(nn_list, label='Neural Network surrogate') + if SUBSETS_ALAMO: + alamo_subset = SUBSETS_ALAMO[i - 1] + alamo_list = alamo(X=alamo_subset[0], P=alamo_subset[1], iters=alamo_subset[2]) + plt.plot(alamo_list, label='ALAMO surrogate') + plt.xlabel('Iteration') + plt.ylabel('Condition number of constraint Jacobian') + plt.title(f'Subset {i}: X={full_subset[0]}, P={full_subset[1]} Pa') + if unsuccessful: + plt.legend(loc="lower right") + else: + plt.legend(loc="upper right") + plt.yscale('log') + plt.grid(True) + plt.tight_layout() + + plt.show() + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Plot subsets") + parser.add_argument("--unsuccessful", action="store_true", help="Whether to use an unsuccessful instance or not") + args = parser.parse_args() + + if args.unsuccessful: + SUBSETS_FULLSPACE = NOSUCCESS_SUBSETS_FULLSPACE + SUBSETS_IMPLICIT = NOSUCCESS_SUBSETS_IMPLICIT + SUBSETS_NN = NOSUCCESS_SUBSETS_NN + fig = plot_subsets(SUBSETS_FULLSPACE, SUBSETS_IMPLICIT, SUBSETS_NN, unsuccessful=True) + else: + SUBSETS_FULLSPACE = SUCCESS_SUBSETS_FULLSPACE + SUBSETS_IMPLICIT = SUCCESS_SUBSETS_IMPLICIT + SUBSETS_ALAMO = SUCCESS_SUBSETS_ALAMO + SUBSETS_NN = SUCCESS_SUBSETS_NN + fig = plot_subsets(SUBSETS_FULLSPACE, SUBSETS_IMPLICIT, SUBSETS_NN, SUBSETS_ALAMO) diff --git a/svi/auto_thermal_reformer/config.py b/svi/auto_thermal_reformer/config.py index f573806..2a89def 100644 --- a/svi/auto_thermal_reformer/config.py +++ b/svi/auto_thermal_reformer/config.py @@ -23,6 +23,83 @@ import itertools import pyomo.environ as pyo from svi.cyipopt import TimedPyomoCyIpoptSolver, Callback +from svi.external import add_external_function_libraries_to_environment +from idaes.core.surrogate.keras_surrogate import KerasSurrogate + + +"""Consistent-signature callbacks to construct each model + +All models returned are initialized and ready to solve. + +Note that these constructors use local imports to avoid circular dependencies. +The right solution here is probably to separate the single-model simulation drivers +(which depend on config) from the model constructing "library" functions (which we +need for these constructor callbacks). + +""" + +def fullspace_constructor(X, P, **kwds): + from svi.auto_thermal_reformer.fullspace_flowsheet import make_optimization_model as create_fullspace_instance + m = create_fullspace_instance(X, P) + return m + + +def alamo_constructor(X, P, **kwds): + from svi.auto_thermal_reformer.alamo_flowsheet import ( + create_instance as create_alamo_instance, + initialize_alamo_atr_flowsheet, + ) + surrogate_fname = kwds.pop("surrogate_fname", None) + m = create_alamo_instance(X, P, surrogate_fname=surrogate_fname) + initialize_alamo_atr_flowsheet(m) + m.fs.reformer_bypass.inlet.temperature.unfix() + m.fs.reformer_bypass.inlet.flow_mol.unfix() + return m + + +def nn_full_constructor(X, P, **kwds): + from svi.auto_thermal_reformer.nn_flowsheet import ( + create_instance as create_nn_instance, + initialize_nn_atr_flowsheet, + ) + surrogate_fname = kwds.pop("surrogate_fname", None) + m = create_nn_instance( + X, + P, + surrogate_fname=surrogate_fname, + formulation=KerasSurrogate.Formulation.FULL_SPACE, + ) + initialize_nn_atr_flowsheet(m) + m.fs.reformer_bypass.inlet.temperature.unfix() + m.fs.reformer_bypass.inlet.flow_mol.unfix() + return m + + +def nn_reduced_constructor(X, P, **kwds): + from svi.auto_thermal_reformer.nn_flowsheet import ( + create_instance as create_nn_instance, + initialize_nn_atr_flowsheet, + ) + surrogate_fname = kwds.pop("surrogate_fname", None) + m = create_nn_instance( + X, + P, + surrogate_fname=surrogate_fname, + formulation=KerasSurrogate.Formulation.REDUCED_SPACE, + ) + initialize_nn_atr_flowsheet(m) + m.fs.reformer_bypass.inlet.temperature.unfix() + m.fs.reformer_bypass.inlet.flow_mol.unfix() + return m + + +def implicit_constructor(X, P, **kwds): + from svi.auto_thermal_reformer.fullspace_flowsheet import make_optimization_model as create_fullspace_instance + from svi.auto_thermal_reformer.implicit_flowsheet import make_implicit + m = create_fullspace_instance(X, P) + add_external_function_libraries_to_environment(m) + m_implicit = make_implicit(m) + return m_implicit filedir = os.path.dirname(__file__) @@ -44,16 +121,25 @@ ] -def get_optimization_solver(options=None): +CONSTRUCTOR_LOOKUP = { + "fullspace": fullspace_constructor, + "implicit": implicit_constructor, + "alamo": alamo_constructor, + "nn-full": nn_full_constructor, + "nn-reduced": nn_reduced_constructor, +} + +def get_optimization_solver(options=None, iters=300, callback=None): # Use cyipopt for everything for Ipopt version consistency among all # formulations #solver = pyo.SolverFactory("cyipopt") # This is a very simple callback we just use to get iteration counts. - cb = Callback() - solver = TimedPyomoCyIpoptSolver(intermediate_callback=cb) + if callback is None: + callback = Callback() + solver = TimedPyomoCyIpoptSolver(intermediate_callback=callback) if options is None: options = {} - solver.config.options["max_iter"] = 300 + solver.config.options["max_iter"] = iters solver.config.options["linear_solver"] = "ma27" solver.config.options["tol"] = 1e-7 solver.config.options["print_user_options"] = "yes" @@ -130,3 +216,14 @@ def get_parameter_samples(args): subset = [int(i) for i in subset] xp_samples = [xp_samples[i] for i in subset] return xp_samples + +def get_plot_argparser(): + argparser = get_argparser() + argparser.add_argument("--show", action="store_true", help="Flag to show the plot") + argparser.add_argument("--no-save", action="store_true", help="Flag to not save the plot") + argparser.add_argument("--plot-fname", default=None, help="Basename for plot file") + argparser.add_argument("--no-legend", action="store_true", help="Flag to exclude a legend") + argparser.add_argument("--title", default=None, help="Plot title") + argparser.add_argument("--show-training-bounds", action="store_true") + argparser.add_argument("--opaque", action="store_true", help="Not transparent") + return argparser diff --git a/svi/auto_thermal_reformer/fullspace_flowsheet.py b/svi/auto_thermal_reformer/fullspace_flowsheet.py index 841cf1a..84c27cd 100644 --- a/svi/auto_thermal_reformer/fullspace_flowsheet.py +++ b/svi/auto_thermal_reformer/fullspace_flowsheet.py @@ -70,6 +70,7 @@ from svi.external import add_external_function_libraries_to_environment from svi.auto_thermal_reformer.reactor_model import add_reactor_model +import svi.auto_thermal_reformer.config as config import argparse @@ -383,6 +384,7 @@ def add_obj_and_constraints(m): ####### OBJECTIVE IS TO MAXIMIZE H2 COMPOSITION IN PRODUCT STREAM ####### m.fs.obj = pyo.Objective( expr=m.fs.product.mole_frac_comp[0, "H2"], sense=pyo.maximize + #expr=-m.fs.product.mole_frac_comp[0, "H2"], sense=pyo.minimize ) # MINIMUM PRODUCT FLOW OF 3500 mol/s IN PRODUCT STREAM @@ -465,10 +467,15 @@ def make_optimization_model(X,P,initialize=True): simulate = not args.optimize if args.optimize: - P = 1600000.0 + P = 1500000.0 X = 0.95 m = make_optimization_model(X, P, initialize=True) - res = pyo.SolverFactory("ipopt").solve(m, tee=True) + m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT) + m.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT) + m.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT) + + solver = config.get_optimization_solver() + res = solver.solve(m, tee=True) elif simulate: P = 1600000.0 diff --git a/svi/auto_thermal_reformer/generate_alamo_surrogate.py b/svi/auto_thermal_reformer/generate_alamo_surrogate.py index 4298c09..c508f49 100644 --- a/svi/auto_thermal_reformer/generate_alamo_surrogate.py +++ b/svi/auto_thermal_reformer/generate_alamo_surrogate.py @@ -124,7 +124,7 @@ def main(): args = argparser.parse_args() - surrogate_fname = os.path.join(args.data_dir, args.surrogate_fname) + surrogate_fname = os.path.join(args.results_dir, args.surrogate_fname) train_plot = os.path.join(args.results_dir, args.train_plot) val_plot = os.path.join(args.results_dir, args.val_plot) diff --git a/svi/auto_thermal_reformer/generate_iterate_data.py b/svi/auto_thermal_reformer/generate_iterate_data.py new file mode 100644 index 0000000..80f5fe2 --- /dev/null +++ b/svi/auto_thermal_reformer/generate_iterate_data.py @@ -0,0 +1,148 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import svi.auto_thermal_reformer.config as config +from svi.cyipopt import FullStateCallback +import pandas as pd +import numpy as np +import pyomo.environ as pyo + + +def main(args): + if args.sample is None: + conversion = 0.94 + pressure = 1550000.0 + else: + xp_samples = config.get_parameter_samples(args) + conversion, pressure = xp_samples[args.sample] + m = config.CONSTRUCTOR_LOOKUP[args.model](conversion, pressure) + m.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT) + m.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT) + + dof_varnames = [ + "fs.reformer_bypass.split_fraction[0.0,bypass_outlet]", + "fs.reformer_mix.steam_inlet_state[0.0].flow_mol", + "fs.feed.properties[0.0].flow_mol", + ] + callback = FullStateCallback( + include_condition=args.include_condition, + include_block_condition=args.include_block_condition, + dof_varnames=dof_varnames, + ) + solver = config.get_optimization_solver(callback=callback) + + solver.config.return_nlp = True + results, nlp = solver.solve(m, tee=True) + + if args.fname is None: + fname = f"{args.model}-iterates" + if args.sample is not None: + fname += f"-{args.sample}" + fname += ".csv" + else: + fname = args.fname + fpath = os.path.join(args.data_dir, fname) + + iterate_data = dict(callback.iterate_data) + # TODO: Should the keys here be updated to indicate that they are primal + # values and residuals? Probably, but I'll deal with that later. + iterate_data.update(callback.primal_values) + iterate_data.update(callback.primal_residuals) + # Make sure nothing pathological happened, like we had a variable named "inf_du" + assert len(iterate_data) == ( + len(callback.iterate_data) + len(callback.primal_values) + len(callback.primal_residuals) + ) + + if args.include_condition: + iterate_data["condition-number"] = list(callback.condition_numbers) + if args.include_block_condition: + iterate_data.update(callback.block_condition_numbers) + + df = pd.DataFrame(iterate_data) + if not args.no_save: + print(f"Saving iterate data to {fpath}") + df.to_csv(fpath) + else: + print("--no-save is set. Not saving iterate data") + + dof_vars = [m.find_component(v) for v in dof_varnames] + from pyomo.common.collections import ComponentSet + dofvarset = ComponentSet(dof_vars) + othervars = [var for var in nlp.get_pyomo_variables() if var not in dofvarset] + var_order = dof_vars + othervars + + from svi.nlp import project_onto, get_reduced_hessian, get_gradient_of_lagrangian + reorder_lbmult = [m.ipopt_zL_out[v] for v in var_order] + reorder_ubmult = [m.ipopt_zU_out[v] for v in var_order] + lbmult = [m.ipopt_zL_out[v] for v in nlp.get_pyomo_variables()] + ubmult = [m.ipopt_zU_out[v] for v in nlp.get_pyomo_variables()] + grad_lag = get_gradient_of_lagrangian(nlp, lbmult, ubmult) + + rh = get_reduced_hessian(nlp, dof_vars, lbmult, ubmult) + eigenvalues, eigenvectors = np.linalg.eig(rh) + print(f"Eigenvalues of reduced Hessian: {eigenvalues}") + sample_suff = "" if args.sample is None else f"-{args.sample}" + rh_fname = f"{args.model}-rh{sample_suff}.npy" + rh_fpath = os.path.join(args.data_dir, rh_fname) + if args.no_save: + print(f"--no-save set. Not saving reduced hessian. Would have saved to {rh_fname}") + else: + np.save(rh_fpath, rh) + #proj_hess = project_onto(rh, [0, 1]) + #eigenvalues, eigenvectors = np.linalg.eig(proj_hess) + + +if __name__ == "__main__": + argparser = config.get_sweep_argparser() + + argparser.add_argument( + "--model", + default="fullspace", + help=( + "Options are 'fullspace', 'implcit', 'alamo', 'nn-full', or 'nn-reduced'." + + " Default is 'fullspace'." + ), + ) + argparser.add_argument( + "--fname", + default=None, + help="Basename for iterate data CSV file. Default is: {model}-iterates-{sample}.csv", + ) + argparser.add_argument( + "--sample", + default=None, + type=int, + help="Index of conversion/pressure parameter sample to use. Default is X=0.95, P=1.55 MPa", + ) + argparser.add_argument( + "--include-condition", + action="store_true", + ) + argparser.add_argument( + "--include-block-condition", + action="store_true", + ) + + args = argparser.parse_args() + if args.subset is not None: + raise RuntimeError("--subset cannot be provided. Use --sample instead") + main(args) diff --git a/svi/auto_thermal_reformer/generate_training_data.py b/svi/auto_thermal_reformer/generate_training_data.py index 10aba3a..8b44d7f 100644 --- a/svi/auto_thermal_reformer/generate_training_data.py +++ b/svi/auto_thermal_reformer/generate_training_data.py @@ -64,8 +64,26 @@ def atr_data_gen( num_samples=600, regular_samples=False, ): - df = {'Fin_CH4':[], 'Tin_CH4':[], 'Fin_H2O':[], 'Conversion': [], 'HeatDuty':[], 'Fout':[], 'Tout':[], 'H2':[], - 'CO':[], 'H2O':[], 'CO2':[], 'CH4':[], 'C2H6':[], 'C3H8':[], 'C4H10':[], 'N2':[], 'O2':[], 'Ar':[]} + df = { + 'Fin_CH4':[], + 'Tin_CH4':[], + 'Fin_H2O':[], + 'Conversion': [], + 'HeatDuty':[], + 'Fout':[], + 'Tout':[], + 'H2':[], + 'CO':[], + 'H2O':[], + 'CO2':[], + 'CH4':[], + 'C2H6':[], + 'C3H8':[], + 'C4H10':[], + 'N2':[], + 'O2':[], + 'Ar':[], + } # Inputs are in the order conversion, F_H2O, F_CH4, T input_ranges = [(0.8, 0.95), (200.0, 350.0), (600.0, 900.0), (600.0, 900.0)] diff --git a/svi/auto_thermal_reformer/get_block_causing_jump.py b/svi/auto_thermal_reformer/get_block_causing_jump.py new file mode 100644 index 0000000..3be4eaa --- /dev/null +++ b/svi/auto_thermal_reformer/get_block_causing_jump.py @@ -0,0 +1,85 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import pyomo.environ as pyo +from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP +from svi.auto_thermal_reformer.fullspace_flowsheet import ( + make_optimization_model, + make_simulation_model, +) + +from pyomo.contrib.incidence_analysis import IncidenceGraphInterface +import svi.auto_thermal_reformer.config as config +import numpy as np + +# SUBSETS_FULLSPACE represents tuples (X,P,iters), where iters = number of iterations +# before the solve experiences a jump in the condition number. +# E.g., Instance X = 0.94, P = 1.94 MPa experiences a jump in iteration 12 + 1 = 13. +# This code exemplifies how to extract the blocks that are causing this jump for this +# specific instance. + +SUBSETS_FULLSPACE = [ + (0.94, 1947379.0, 12), +] + +def calculate_condition_number(m, condition_number_threshold = 1e6): + result = {"block": [], "condition number": []} + nlp = PyomoNLP(m) + igraph = IncidenceGraphInterface(m, include_inequality=False) + vblocks, cblocks = igraph.block_triangularize() + for i, (vblock, cblock) in enumerate(zip(vblocks, cblocks)): + submatrix = nlp.extract_submatrix_jacobian(vblock, cblock) + cond = np.linalg.cond(submatrix.toarray()) + if cond > condition_number_threshold: + result["block"].append(i) + result["condition number"].append(cond) + return result + +def full_space(X=0.94, P=1947379.0, iters=300): + m = make_optimization_model(X, P) + solver = config.get_optimization_solver(iters=iters) + print(f"Solving sample with X={X}, P={P}") + results = solver.solve(m, tee=True) + m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].fix(m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].value) + m.fs.reformer_mix.steam_inlet.flow_mol.fix(m.fs.reformer_mix.steam_inlet.flow_mol[0].value) + m.fs.feed.outlet.flow_mol.fix(m.fs.feed.outlet.flow_mol[0].value) + result = calculate_condition_number(m) + return result + +def id_blocks_high_condition_number(dict1, dict2): + increased_blocks = [] + for block1, cond1, cond2 in zip(dict1['block'], dict1['condition number'], dict2['condition number']): + if cond2 > cond1: + increased_blocks.append(block1) + return increased_blocks + +def main(): + for i, (full_subset) in enumerate(zip(SUBSETS_FULLSPACE), 1): + before_jump = full_space(X=full_subset[0][0], P=full_subset[0][1], iters=full_subset[0][2]) + after_jump = full_space(X=full_subset[0][0], P=full_subset[0][1], iters=full_subset[0][2]+1) + block_to_analyze = id_blocks_high_condition_number(before_jump, after_jump) + # Maybe also print equations and variables present in that block? For now, all we need is the number. + # With this number we will proceed to identify and print variable values participating in this block. + print(f"Block(s) causing the instance to enter a region of high condition number: {block_to_analyze}.") + +if __name__ == "__main__": + main() diff --git a/svi/auto_thermal_reformer/identify_vars.py b/svi/auto_thermal_reformer/identify_vars.py new file mode 100644 index 0000000..1f8c822 --- /dev/null +++ b/svi/auto_thermal_reformer/identify_vars.py @@ -0,0 +1,143 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import pyomo.environ as pyo +from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP +from svi.auto_thermal_reformer.fullspace_flowsheet import ( + make_optimization_model, + make_simulation_model, +) +import pandas as pd +from svi.auto_thermal_reformer.implicit_flowsheet import make_implicit +from svi.external import add_external_function_libraries_to_environment +from pyomo.contrib.pynumero.interfaces.external_pyomo_model import ExternalPyomoModel +from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock +from pyomo.contrib.pynumero.algorithms.solvers.implicit_functions import ( + CyIpoptSolverWrapper +) + +from pyomo.core.expr.visitor import identify_variables +from pyomo.util.subsystems import create_subsystem_block +from pyomo.contrib.incidence_analysis import IncidenceGraphInterface +import svi.auto_thermal_reformer.config as config +import numpy as np +import argparse + +# SUBSETS_FULLSPACE represents tuples (X,P,iters), where iters = number of iterations +# before the solve experiences a jump in the condition number. +# E.g., Instance X = 0.94, P = 1.94 MPa experiences a jump in iteration 12 + 1 = 13. + +SUBSETS_FULLSPACE = [ + (0.94, 1947379.0, 12), +] + +# We want to see what variable values are causing this jump in condition number. +# For this specific subset, jump in condition number occurs only in block 524. + +def get_vars_related_to_blocks(m, blocks=[524]): + input_and_block_vars = [] + nlp = PyomoNLP(m) + igraph = IncidenceGraphInterface(m, include_inequality=False) + vblocks, cblocks = igraph.block_triangularize() + for block in blocks: + block_vars = {"Block": block, "Variable": [], "Value": []} + vblock = vblocks[block] + cblock = cblocks[block] + + for con in cblock: + for var in identify_variables(con.expr): + if 'params' not in var.name: + block_vars["Variable"].append(var.name) + block_vars["Value"].append(var.value) + + input_and_block_vars.append(block_vars) + return input_and_block_vars + +def full_space(X=0.94, P=1947379.0, iters=300, blocks=[524]): + m = make_optimization_model(X, P) + solver = config.get_optimization_solver(iters=iters) + print(f"Solving sample with X={X}, P={P}") + results = solver.solve(m, tee=True) + m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].fix(m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].value) + m.fs.reformer_mix.steam_inlet.flow_mol.fix(m.fs.reformer_mix.steam_inlet.flow_mol[0].value) + m.fs.feed.outlet.flow_mol.fix(m.fs.feed.outlet.flow_mol[0].value) + dicts = get_vars_related_to_blocks(m, blocks = blocks) + return dicts + +def implicit(X=0.94, P=1947379.0, iters=300, blocks = [524]): + m = make_optimization_model(X, P) + add_external_function_libraries_to_environment(m) + m_implicit = make_implicit(m) + solver = config.get_optimization_solver(iters=iters) + print(f"Solving sample with X={X}, P={P}") + results = solver.solve(m_implicit, tee=True) + m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].fix(m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].value) + m.fs.reformer_mix.steam_inlet.flow_mol.fix(m.fs.reformer_mix.steam_inlet.flow_mol[0].value) + m.fs.feed.outlet.flow_mol.fix(m.fs.feed.outlet.flow_mol[0].value) + dicts = get_vars_related_to_blocks(m, blocks = blocks) + return dicts + +# The function below is written under the premise that a jump in condition number is due to +# a significant jump in some variables values. The magnitude of this jump for each variable +# is unknown; percentage_diff argument can me modified according to the user's experience. + +def variable_with_jump(dict1, dict2, percentage_diff = 10): + variables_with_jump = {"Variable":[], "Value_initial":[], "Value_final":[]} + for var, value1, value2 in zip(dict1["Variable"], dict1["Value"], dict2["Value"]): + if value2 != value1: + percentage_difference = abs((value2 - value1) / value1) * 100 + if percentage_difference > percentage_diff: + variables_with_jump["Variable"].append(var) + variables_with_jump["Value_initial"].append(value1) + variables_with_jump["Value_final"].append(value2) + return variables_with_jump + +def dict_to_csv(data, filename): + df = pd.DataFrame.from_dict(data) + df.to_csv(filename, index=False) + +def main(blocks, percentage_diff, use_full_space, full_subsets): + for full_subset in full_subsets: + if use_full_space: + before_jump = full_space(X=full_subset[0], P=full_subset[1], iters=full_subset[2], blocks = blocks) + after_jump = full_space(X=full_subset[0], P=full_subset[1], iters=full_subset[2] + 1, blocks = blocks) + else: + before_jump = implicit(X=full_subset[0], P=full_subset[1], iters=full_subset[2], blocks = blocks) + after_jump = implicit(X=full_subset[0], P=full_subset[1], iters=full_subset[2] + 1, blocks = blocks) + + for before, after in zip(before_jump, after_jump): + block_number = before["Block"] + variables_with_jump = variable_with_jump(before, after, percentage_diff) + filename = f'vars_jump_more_than_{percentage_diff}pt_block_{block_number}.csv' + dict_to_csv(variables_with_jump, filename) + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--blocks", nargs="+", type=int, default=[524], help="List of blocks") + parser.add_argument("--percentage_diff", type=float, default=10, help="Threshold for percentage difference") + parser.add_argument("--use_full_space", action="store_true", help="Use full_space function") + args = parser.parse_args() + + full_subsets = SUBSETS_FULLSPACE + main(args.blocks, args.percentage_diff, args.use_full_space, full_subsets) + + diff --git a/svi/auto_thermal_reformer/implicit_flowsheet.py b/svi/auto_thermal_reformer/implicit_flowsheet.py index 50316e4..d6f1677 100644 --- a/svi/auto_thermal_reformer/implicit_flowsheet.py +++ b/svi/auto_thermal_reformer/implicit_flowsheet.py @@ -38,6 +38,7 @@ from idaes.models.properties.modular_properties import GenericParameterBlock from idaes.core.util.model_statistics import degrees_of_freedom +import idaes from idaes.models.unit_models import ( Mixer, Heater, @@ -118,8 +119,12 @@ def make_implicit(m, **kwds): # coefficients with values of zero were handled. This should be fixed in # recent Pyomo main, but we leave this explicit filtering in here as we # often run using the latest release. - to_exclude.add(m.fs.reformer.lagrange_mult[0, "N"]) - to_exclude.add(m.fs.reformer.lagrange_mult[0, "Ar"]) + # + # NOTE: As of IDAES 2.5, lagrange multipliers for inert components are not + # included in the Gibbs reactor. + if (idaes.ver.package_version.major, idaes.ver.package_version.minor) < (2, 5): + to_exclude.add(m.fs.reformer.lagrange_mult[0, "N"]) + to_exclude.add(m.fs.reformer.lagrange_mult[0, "Ar"]) external_vars = [var for var in reformer_igraph.variables if var not in to_exclude] external_var_set = ComponentSet(external_vars) diff --git a/svi/auto_thermal_reformer/makefile b/svi/auto_thermal_reformer/makefile index 8a38ef3..45555bd 100644 --- a/svi/auto_thermal_reformer/makefile +++ b/svi/auto_thermal_reformer/makefile @@ -3,7 +3,7 @@ SUBSET = "0,2,4,6,7,8,14,18,21,24,25,26,27,30,32,35,36,37,38,39,41,42,44,45,47" train: python generate_training_data.py --regular-samples python generate_alamo_surrogate.py data/training-data.csv - python nn_tuning_training.py data/training-data.csv + python nn_tuning_training.py data/training-data.csv --tune sweep: python run_fullspace_sweep.py @@ -62,3 +62,33 @@ plot: --validation-fpath=data/implicit-sweep-validation.csv \ --feastol=1e-5 \ --title="Implicit function" + +plot-error: + python plot_surrogate_error.py \ + data/alamo-sweep-surrogate-error.csv \ + --no-legend \ + --opaque \ + --title="ALAMO" + python plot_surrogate_error.py \ + data/nn-sweep-full-surrogate-error.csv \ + --opaque \ + --title="Neural network" + +plot-error-png: + python plot_surrogate_error.py \ + data/alamo-sweep-surrogate-error.csv \ + --no-legend \ + --opaque \ + --title="ALAMO" \ + --plot-fname=alamo-sweep-surrogate-error.png + python plot_surrogate_error.py \ + data/nn-sweep-full-surrogate-error.csv \ + --opaque \ + --title="Neural network" \ + --plot-fname=nn-sweep-full-surrogate-error.png + +plot-nn-state15: + python plot_state.py data-202405/nn-reduced-iterates-15.csv,data-202405/nn-full-iterates-15.csv split-fraction --state2=methane-inlet --rh=data-202405/nn-full-rh-15.npy --output-fname=nn-states-15.pdf + +plot-nn-state14: + python plot_state.py data-202405/nn-reduced-iterates-14.csv,data-202405/nn-full-iterates-14.csv split-fraction --state2=methane-inlet --rh=data-202405/nn-full-rh-14.npy --output-fname=nn-states-14.pdf diff --git a/svi/auto_thermal_reformer/nn_tuning_training.py b/svi/auto_thermal_reformer/nn_tuning_training.py index 4cfccb0..d92edf7 100644 --- a/svi/auto_thermal_reformer/nn_tuning_training.py +++ b/svi/auto_thermal_reformer/nn_tuning_training.py @@ -136,6 +136,14 @@ def gibbs_to_nn( keras_surrogate.save_to_folder(surrogate_fname) print(f"Saved NN surrogate model to {surrogate_fname}") + print( + f"Parameters of model that was just saved:" + f"activation={activation}," + f"optimizer={optimizer}," + f"n_hidden_layers={n_hidden_layers}," + f"n_nodes_per_payer={n_nodes_per_layer}," + ) + t1 = time.time() total_time = t1 - t0 print("Total time: ", total_time) diff --git a/svi/auto_thermal_reformer/plot_condition_num_grid.py b/svi/auto_thermal_reformer/plot_condition_num_grid.py new file mode 100644 index 0000000..335b2cd --- /dev/null +++ b/svi/auto_thermal_reformer/plot_condition_num_grid.py @@ -0,0 +1,118 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# __________________________________________________________________________ + +import os +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.patches import Patch +import seaborn as sns +import svi.auto_thermal_reformer.config as config +from matplotlib.colors import LogNorm + +def plot_cond_num( + df, + legend=True, + show_training_bounds=False, +): + # TODO: Make this (as well as font.family) a global option set + # when we get the argparser? Or maybe just make this configurable from the + # command line with a global default? + # But the same font size might not be appropriate for all plots. + plt.rcParams["font.size"] = 20 + fig = plt.figure() + + conversion_list = list(sorted(set(df["X"]))) + pressure_list = list(sorted(set(df["P"]))) + n_conversion = len(conversion_list) + n_pressure = len(pressure_list) + cond_num_lookup = { + (df["X"][i], df["P"][i]): df["Condition Number"][i] + for i in range(len(df)) + } + cond_num_array = [ + [cond_num_lookup[conversion_list[j], pressure_list[i]] for i in range(n_pressure)] + for j in range(n_conversion) + ] + + cbar = legend + cbar_kws = dict(label="Condition Number") + ax = sns.heatmap( + cond_num_array, + cbar=cbar, + norm=LogNorm(vmin=1e19, vmax=1e27), + cbar_kws=cbar_kws, + square=True, + cmap="Reds", + ) + + xtick_positions = [i+0.5 for i in range(n_pressure)] + ytick_positions = [i+0.5 for i in range(n_conversion)] + xtick_labels = ["%1.2f" % (p / 1e6) if i%2 else "" for i, p in enumerate(pressure_list)] + ytick_labels = ["%1.2f" % x if i%2 else "" for i, x in enumerate(conversion_list)] + ax.set_xticks(xtick_positions, labels=xtick_labels) + ax.set_yticks(ytick_positions, labels=ytick_labels, rotation=0) + ax.set_xlabel("Pressure (MPa)") + ax.set_ylabel("Conversion") + + plt.gca().invert_yaxis() + ax.set_facecolor("black") + + return fig, ax + + +def main(args): + df = pd.read_csv(args.cond_num_fpath) + + fig, ax = plot_cond_num( + df, + legend=not args.no_legend, + show_training_bounds=args.show_training_bounds, + ) + + if args.title is not None: + ax.set_title(args.title) + + fig.tight_layout() + + if not args.no_save: + if args.plot_fname is None: + plot_fname = os.path.basename(args.cond_num_fpath) + data_ext = "." + plot_fname.split(".")[-1] + ext_len = len(data_ext) + plot_fname = plot_fname[:-ext_len] + "-condition.pdf" + else: + plot_fname = args.plot_fname + + plot_fpath = os.path.join(args.results_dir, plot_fname) + print(f"Saving figure to {plot_fpath}") + fig.savefig(plot_fpath, transparent=not args.opaque) + + if args.show: + plt.show() + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + argparser.add_argument( + "cond_num_fpath", help="Path to CSV file containing condition numbers to plot" + ) + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/plot_implicit_imat.py b/svi/auto_thermal_reformer/plot_implicit_imat.py index 3c2e50b..ebbd49b 100644 --- a/svi/auto_thermal_reformer/plot_implicit_imat.py +++ b/svi/auto_thermal_reformer/plot_implicit_imat.py @@ -37,6 +37,7 @@ from idaes.models.properties.modular_properties import GenericParameterBlock from idaes.core.util.model_statistics import degrees_of_freedom +import idaes from idaes.models.unit_models import ( Mixer, Heater, @@ -112,8 +113,12 @@ def make_implicit(m): # coefficients with values of zero were handled. This should be fixed in # recent Pyomo main, but we leave this explicit filtering in here as we # often run using the latest release. - to_exclude.add(m.fs.reformer.lagrange_mult[0, "N"]) - to_exclude.add(m.fs.reformer.lagrange_mult[0, "Ar"]) + # + # NOTE: As of IDAES v... (2.5 or so), these inert components don't exist in + # the model. + if (idaes.ver.package_version.major, idaes.ver.package_version.minor) < (2, 5): + to_exclude.add(m.fs.reformer.lagrange_mult[0, "N"]) + to_exclude.add(m.fs.reformer.lagrange_mult[0, "Ar"]) external_vars = [var for var in reformer_igraph.variables if var not in to_exclude] external_var_set = ComponentSet(external_vars) diff --git a/svi/auto_thermal_reformer/plot_iterate_block_condition_number.py b/svi/auto_thermal_reformer/plot_iterate_block_condition_number.py new file mode 100644 index 0000000..c7f1cd1 --- /dev/null +++ b/svi/auto_thermal_reformer/plot_iterate_block_condition_number.py @@ -0,0 +1,118 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import svi.auto_thermal_reformer.config as config +import pandas as pd +import matplotlib.pyplot as plt + + +LABEL_LOOKUP = { + "inf_pr": "Primal infeasibility", + "inf_du": "Dual infeasibility", + "condition-number": "Condition number", + "fullspace": "Full-space", + "implicit": "Implicit function", + "nn-full": "Neural network", + "alamo": "ALAMO surrogate", +} + + +FONTSIZE = 16 + + +def plot_trajectory( + df, + keys, + labels=None, + fig=None, + ax=None, + iter_lim=None, + legend=True, +): + plt.rcParams["font.size"] = FONTSIZE + if fig is None or ax is None: + fig, ax = plt.subplots() + + iterations = list(range(len(df))) + + for i, key in enumerate(keys): + if labels is None: + label = LABEL_LOOKUP[key] if key in LABEL_LOOKUP else key + else: + label = labels[i] + ax.plot( + iterations, + list(df[key]), + label=label, + linewidth=2, + ) + if legend: + ax.legend() + ax.set_yscale("log") + ax.xaxis.set_tick_params(length=0) + ax.yaxis.set_tick_params(length=0) + ax.set_xlabel("Iteration number") + if iter_lim is not None: + ax.set_xlim(None, iter_lim) + return fig, ax + + +def main(args): + fpath = args.fpath + df = pd.read_csv(fpath) + keys = [colname for colname in df.columns if colname.startswith("block") and colname.endswith("cond")] + + plot_trajectory( + df, + keys, + iter_lim=args.iter_lim, + legend=not args.no_legend, + ) + + if args.show: + plt.show() + + if not args.no_save: + fig.tight_layout() + # Assume file name is NAME.ext + name = os.path.basename(args.fpath).split(".")[0] + fname = name + "-block-condition.pdf" + fpath = os.path.join(args.results_dir, fname) + fig.savefig(fpath, transparent=not args.opaque) + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + + argparser.add_argument( + "fpath", + help="File containing iterate data to plot", + ) + argparser.add_argument( + "--iter-lim", + type=int, + default=None, + help="Upper bound on x-axis", + ) + + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/plot_iterate_condition_number.py b/svi/auto_thermal_reformer/plot_iterate_condition_number.py new file mode 100644 index 0000000..afe74d6 --- /dev/null +++ b/svi/auto_thermal_reformer/plot_iterate_condition_number.py @@ -0,0 +1,143 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import svi.auto_thermal_reformer.config as config +import pandas as pd +import matplotlib.pyplot as plt + + +LABEL_LOOKUP = { + "inf_pr": "Primal infeasibility", + "inf_du": "Dual infeasibility", + "condition-number": "Condition number", + "fullspace": "Full-space", + "implicit": "Implicit function", + "nn-full": "Neural network", + "alamo": "ALAMO surrogate", +} + + +FONTSIZE = 16 + + +def plot_trajectory( + df, + keys, + labels=None, + fig=None, + ax=None, + iter_lim=None, +): + plt.rcParams["font.size"] = FONTSIZE + if fig is None or ax is None: + fig, ax = plt.subplots() + + iterations = list(range(len(df))) + + for i, key in enumerate(keys): + if labels is None: + label = LABEL_LOOKUP[key] if key in LABEL_LOOKUP else key + else: + label = labels[i] + ax.plot( + iterations, + list(df[key]), + label=label, + linewidth=2, + ) + if not args.no_legend: + ax.legend() + ax.set_yscale("log") + ax.xaxis.set_tick_params(length=0) + ax.yaxis.set_tick_params(length=0) + ax.set_xlabel("Iteration number") + if iter_lim is not None: + ax.set_xlim(None, iter_lim) + return fig, ax + + +def main(args): + fpaths = args.fpaths.split(",") + models = [] + for fpath in fpaths: + model = None + for name in config.CONSTRUCTOR_LOOKUP: + if name in os.path.basename(fpath): + model = name + break + if model is None: + raise RuntimeError( + f"Cannot infer model from fpath {fpath}" + ) + models.append(model) + + dfs = [pd.read_csv(fpath) for fpath in fpaths] + #keys = [key for key in args.keys.split(",") if key != ""] + keys = ["condition-number"] + labels = [LABEL_LOOKUP[model] for model in models] + + plt.rcParams["font.size"] = FONTSIZE + fig, ax = plt.subplots() + for i, df in enumerate(dfs): + plot_trajectory( + df, + keys, + # Use the labels arg here to override the default, which is to lookup + # labels by keys. + labels=[labels[i]], + fig=fig, + ax=ax, + iter_lim=args.iter_lim, + ) + + if args.show: + plt.show() + + if not args.no_save: + fig.tight_layout() + # Assume file name is NAME.ext + if len(models) == 1: + name = os.path.basename(args.fpath).split(".")[0] + fname = name + "-condition.pdf" + fpath = os.path.join(args.results_dir, fname) + else: + fname = "-".join(models) + "-iterates-condition.pdf" + fpath = os.path.join(args.results_dir, fname) + fig.savefig(fpath, transparent=not args.opaque) + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + + argparser.add_argument( + "fpaths", + help="Comma-separated list of files containing iterate data to plot", + ) + argparser.add_argument( + "--iter-lim", + type=int, + default=None, + help="Upper bound on x-axis", + ) + + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/plot_primal_dual_infeas.py b/svi/auto_thermal_reformer/plot_primal_dual_infeas.py new file mode 100644 index 0000000..4172458 --- /dev/null +++ b/svi/auto_thermal_reformer/plot_primal_dual_infeas.py @@ -0,0 +1,84 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import svi.auto_thermal_reformer.config as config +import pandas as pd +import matplotlib.pyplot as plt + + +LABEL_LOOKUP = { + "inf_pr": "Primal infeasibility", + "inf_du": "Dual infeasibility", +} + + +def plot_trajectory( + df, + keys, + labels=None, +): + plt.rcParams["font.size"] = 16 + fig, ax = plt.subplots() + + iterations = list(range(len(df))) + + for i, key in enumerate(keys): + label = LABEL_LOOKUP[key] if key in LABEL_LOOKUP else key + ax.plot( + iterations, + list(df[key]), + label=label, + linewidth=2, + ) + ax.legend() + ax.set_yscale("log") + ax.xaxis.set_tick_params(length=0) + ax.yaxis.set_tick_params(length=0) + ax.set_xlabel("Iteration number") + return fig, ax + + +def main(args): + df = pd.read_csv(args.fpath) + #keys = [key for key in args.keys.split(",") if key != ""] + keys = ["inf_pr", "inf_du"] + fig, ax = plot_trajectory(df, keys) + + if args.show: + plt.show() + + if not args.no_save: + fig.tight_layout() + # Assume file name is NAME.ext + name = os.path.basename(args.fpath).split(".")[0] + fname = name + "-pdinfeas.pdf" + fpath = os.path.join(args.results_dir, fname) + fig.savefig(fpath, transparent=not args.opaque) + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + + argparser.add_argument("fpath", help="File containing iterate data to plot") + + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/plot_residual_histogram.py b/svi/auto_thermal_reformer/plot_residual_histogram.py new file mode 100644 index 0000000..99d1c9a --- /dev/null +++ b/svi/auto_thermal_reformer/plot_residual_histogram.py @@ -0,0 +1,92 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import svi.auto_thermal_reformer.config as config +import pandas as pd +import matplotlib.pyplot as plt +import numpy as np +from pyomo.core.base.constraint import Constraint + + +def main(args): + df = pd.read_csv(args.fpath) + + model = None + for key in config.CONSTRUCTOR_LOOKUP: + if key in os.path.basename(args.fpath): + model = key + if model is None: + raise RuntimeError("Could not infer model from filename") + # Just choose some dummy values of conversion and pressure. We won't + # actually use this model numerically + X = 0.94 + P = 1550000.0 + m = config.CONSTRUCTOR_LOOKUP[model](X, P) + connames = [con.name for con in m.component_data_objects(Constraint)] + + con_infeas_count = {} + for name in connames: + resid_array = np.array(df[name]) + violated = resid_array > args.infeas_threshold + n_violated = sum(violated) + con_infeas_count[name] = n_violated + + constraints_by_infeas = sorted( + connames, + reverse=True, + key=lambda name: con_infeas_count[name], + ) + infeas_counts = [con_infeas_count[name] for name in constraints_by_infeas] + constraint_indices = list(range(len(constraints_by_infeas))) + + plt.rcParams["font.size"] = 20 + fig, ax = plt.subplots() + ax.bar(constraint_indices, infeas_counts) + ax.set_ylabel("Number of\niterations", rotation=0) + ax.set_xlabel("Constraint indices") + ax.yaxis.set_label_coords(-0.3, 0.5) + ax.xaxis.set_tick_params(length=0) + ax.yaxis.set_tick_params(length=0) + + w, h = fig.get_size_inches() + fig.set_size_inches(w*1.3, h) + + fig.tight_layout() + if args.show: + plt.show() + + if not args.no_save: + # Assume file name is NAME.ext + name = os.path.basename(args.fpath).split(".")[0] + fname = name + "-resid-histogram.pdf" + fpath = os.path.join(args.results_dir, fname) + fig.savefig(fpath, transparent=not args.opaque) + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + + argparser.add_argument("fpath", help="File containing iterate data to plot") + argparser.add_argument("--infeas-threshold", type=float, default=1e-3) + + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/plot_state.py b/svi/auto_thermal_reformer/plot_state.py new file mode 100644 index 0000000..2af2fc0 --- /dev/null +++ b/svi/auto_thermal_reformer/plot_state.py @@ -0,0 +1,287 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import svi.auto_thermal_reformer.config as config +from svi.nlp import project_onto +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.patches import Ellipse +import numpy as np + + +"""Script for plotting comparisons of state trajectories between different methods +""" + + +LABEL_LOOKUP = { + "inf_pr": "Primal infeasibility", + "inf_du": "Dual infeasibility", + "fs.reformer_bypass.split_fraction[0.0,bypass_outlet]": "Split fraction", + "fs.reformer_mix.steam_inlet_state[0.0].flow_mol": "Steam inlet", + "fs.feed.properties[0.0].flow_mol": "Natural gas inlet", +} + +KEY_LOOKUP = { + # These are names that we can use in the CLI instead of writing out the keys + "split-fraction": "fs.reformer_bypass.split_fraction[0.0,bypass_outlet]", + "steam-inlet": "fs.reformer_mix.steam_inlet_state[0.0].flow_mol", + "methane-inlet": "fs.feed.properties[0.0].flow_mol", +} + +INDEX_LOOKUP = { + "fs.reformer_bypass.split_fraction[0.0,bypass_outlet]": 0, + "fs.reformer_mix.steam_inlet_state[0.0].flow_mol": 1, + "fs.feed.properties[0.0].flow_mol": 2, +} + + +plt.rcParams["font.size"] = 12 +plt.rcParams["font.family"] = "serif" + + +def get_label(fpath): + fname = os.path.basename(fpath) + if "nn-full" in fname: + return "Full-space NN" + elif "nn-reduced" in fname: + return "Reduced-space NN" + elif "fullspace" in fname: + return "Full-space" + elif "implicit" in fname: + return "Implicit" + elif "alamo" in fname: + return "ALAMO" + else: + raise NotImplementedError(f"Filepath {fpath} could not be recognized") + + +def plot_trajectory( + df, + keys, + # Overrides default labels so we can label trajectories by file rather than by key + labels=None, + fig_ax=None, + logscale=False, +): + fig, ax = plt.subplots() if fig_ax is None else fig_ax + + iterations = list(range(len(df))) + + if labels is None: + labels = [LABEL_LOOKUP.get(key, key) for key in keys] + + for i, key in enumerate(keys): + label = labels[i] + ax.plot( + iterations, + list(df[key]), + label=label, + linewidth=2, + ) + ax.legend() + if logscale: + ax.set_yscale("log") + ax.xaxis.set_tick_params(length=0) + ax.yaxis.set_tick_params(length=0) + ax.set_xlabel("Iteration number") + return fig, ax + + +def plot_2d_state( + df, + state1, + state2, + fig_ax=None, + label=None, + state1label=None, + state2label=None, + # Include initial and final points + include_points=True, +): + fig, ax = plt.subplots() if fig_ax is None else fig_ax + state1label = LABEL_LOOKUP.get(state1, state1) if state1label is None else state1label + state2label = LABEL_LOOKUP.get(state2, state2) if state2label is None else state2label + + xdata = list(df[state1]) + ydata = list(df[state2]) + ax.plot( + xdata, + ydata, + label=label, + ) + ax.set_xlabel(state1label) + ax.set_ylabel(state2label) + + ax.legend() + return fig, ax + + +def get_opt_xy(df, state1, state2): + xdata = list(df[state1]) + ydata = list(df[state2]) + finalx = xdata[-1] + finaly = ydata[-1] + return (finalx, finaly) + + +def plot_points( + df, + state1, + state2, + fig_ax=None, + include_final=True, +): + fig, ax = plt.subplots() if fig_ax is None else fig_ax + xdata = list(df[state1]) + ydata = list(df[state2]) + initx = xdata[0] + inity = ydata[0] + finalx = xdata[-1] + finaly = ydata[-1] + ax.scatter( + [initx], + [inity], + marker='.', + s=100, + label="Initial", + color="black", + # Sufficiently high zorder that we put points on top of lines + zorder=10, + ) + ax.scatter( + [finalx], + [finaly], + marker='*', + s=100, + label="Optimal", + color="black", + zorder=10, + ) + ax.legend() + return fig, ax + + +def plot_contours( + rh, + center=(0.0, 0.0), + fig_ax=None, + levels=None, +): + fig, ax = plt.subplots() if fig_ax is None else fig_ax + #levels = [1.0, 2.0, 3.0] if levels is None else levels + levels = [1.0] if levels is None else levels + evals, evecs = np.linalg.eig(rh) + # TODO: Make sure RH is 2x2? + if not np.all(evals > 0): + raise ValueError("Reduced hessian is not positive definite") + for l in levels: + width = 2 * (l / evals[0])**0.5 + height = 2 * (l / evals[1])**0.5 + angle = np.rad2deg(np.arctan2(*evecs[1])) + ellipse = Ellipse( + center, + width=width, + height=height, + angle=angle, + edgecolor=tuple([0.7]*3), + facecolor="none", + ) + ax.add_patch(ellipse) + + +def main(args): + fpaths = args.fpaths.split(",") + dfs = [pd.read_csv(fpath) for fpath in fpaths] + + # TODO: Handle multiple keys in a trajectory plot? + keys = [args.state] + # Get keys if any shorthands were used + keys = [KEY_LOOKUP.get(key, key) for key in keys] + fig, ax = plt.subplots() + if args.state2 is None: + for i, df in enumerate(dfs): + label = get_label(fpaths[i]) + # I know that we're only plotting one key + labels = [label] + plot_trajectory(df, keys, labels=labels, fig_ax=(fig, ax)) + output_fname = "state-trajectory.pdf" if args.output_fname is None else args.output_fname + else: + for i, df in enumerate(dfs): + label = get_label(fpaths[i]) + state2 = KEY_LOOKUP.get(args.state2, args.state2) + plot_2d_state( + df, + keys[0], + state2, + fig_ax=(fig, ax), + label=label, + ) + # We plot the successful trajectory second + plot_points( + dfs[1], + keys[0], + state2, + fig_ax=(fig, ax), + ) + output_fname = "state.pdf" if args.output_fname is None else args.output_fname + if args.rh is not None: + coords = (INDEX_LOOKUP[keys[0]], INDEX_LOOKUP[state2]) + rh = np.load(args.rh) + proj_rh = project_onto(rh, coords) + levels = [0.00625, 0.0125, 0.025, 0.05, 0.1, 0.2] + center = get_opt_xy(dfs[1], keys[0], state2) + plot_contours( + proj_rh, + center=center, + fig_ax=(fig, ax), + levels=levels, + ) + #ax.set_xlim(0.0, 0.6) + #ax.set_ylim(1000, 1400) + ax.legend(loc="upper left") + + fig.tight_layout() + + if args.show: + plt.show() + + if not args.no_save: + # Assume file name is NAME.ext + output_fpath = os.path.join(args.results_dir, output_fname) + fig.savefig(output_fpath, transparent=not args.opaque) + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + + # TODO: Get fpaths + argparser.add_argument( + "fpaths", + help="Comma-separated list of files containing states to plot", + ) + argparser.add_argument("state", help="State to plot. This can be anything we track, e.g. 'obj_value'") + argparser.add_argument("--state2", help="Optional second state to plot", default=None) + argparser.add_argument("--rh", default=None, help=".npy file containing reduced hessian") + argparser.add_argument("--output-fname", default=None) + + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/plot_sum_constraint_residual.py b/svi/auto_thermal_reformer/plot_sum_constraint_residual.py new file mode 100644 index 0000000..81cd18e --- /dev/null +++ b/svi/auto_thermal_reformer/plot_sum_constraint_residual.py @@ -0,0 +1,107 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# ___________________________________________________________________________ + +import os +import svi.auto_thermal_reformer.config as config +import pandas as pd +import matplotlib.pyplot as plt + + +LABEL_LOOKUP = { + "inf_pr": "Primal infeasibility", + "inf_du": "Dual infeasibility", + "fs.reformer.control_volume.properties_out[0.0].sum_mole_frac_out": "Reformer sum mole fraction", + "fs.reformer_recuperator.hot_side.properties_out[0.0].sum_mole_frac_out": "Recuperator sum mole fraction", +} + + +def plot_trajectory( + df, + keys, + labels=None, + iter_lim=None, + legend=True, +): + plt.rcParams["font.size"] = 16 + fig, ax = plt.subplots() + + iterations = list(range(len(df))) + + w, h = fig.get_size_inches() + fig.set_size_inches(w*1.5, h) + for i, key in enumerate(keys): + label = LABEL_LOOKUP[key] if key in LABEL_LOOKUP else key + ax.plot( + iterations, + list(df[key]), + label=label, + linewidth=2, + ) + if legend: + ax.legend(loc=(0.6, 0.5)) + ax.set_yscale("log") + ax.xaxis.set_tick_params(length=0) + ax.yaxis.set_tick_params(length=0) + ax.set_xlabel("Iteration number") + if iter_lim is not None: + ax.set_xlim(-1, iter_lim) + return fig, ax + + +def main(args): + df = pd.read_csv(args.fpath) + #keys = [key for key in args.keys.split(",") if key != ""] + keys = [ + "fs.reformer.control_volume.properties_out[0.0].sum_mole_frac_out", + "fs.reformer_recuperator.hot_side.properties_out[0.0].sum_mole_frac_out", + ] + fig, ax = plot_trajectory( + df, + keys, + iter_lim=args.iter_lim, + legend=not args.no_legend, + ) + + fig.tight_layout() + if args.show: + plt.show() + + if not args.no_save: + # Assume file name is NAME.ext + name = os.path.basename(args.fpath).split(".")[0] + fname = name + "-sum-mole-frac.pdf" + fpath = os.path.join(args.results_dir, fname) + fig.savefig(fpath, transparent=not args.opaque) + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + + argparser.add_argument("fpath", help="File containing iterate data to plot") + argparser.add_argument( + "--iter-lim", + type=int, + default=None, + help="Upper bound on x-axis", + ) + + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/plot_surrogate_error.py b/svi/auto_thermal_reformer/plot_surrogate_error.py new file mode 100644 index 0000000..f3e5ea8 --- /dev/null +++ b/svi/auto_thermal_reformer/plot_surrogate_error.py @@ -0,0 +1,136 @@ +# ___________________________________________________________________________ +# +# Surrogate vs. Implicit: Experiments comparing nonlinear optimization +# formulations +# +# Copyright (c) 2023. Triad National Security, LLC. All rights reserved. +# +# This program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National Nuclear +# Security Administration. All rights in the program are reserved by Triad +# National Security, LLC, and the U.S. Department of Energy/National Nuclear +# Security Administration. The Government is granted for itself and others +# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +# in this material to reproduce, prepare derivative works, distribute copies +# to the public, perform publicly and display publicly, and to permit others +# to do so. +# +# This software is distributed under the 3-clause BSD license. +# __________________________________________________________________________ + +import os +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.patches import Patch +import seaborn as sns +import svi.auto_thermal_reformer.config as config + + +def plot_surrogate_error( + df, + legend=True, + show_training_bounds=False, +): + # TODO: Make this (as well as font.family) a global option set + # when we get the argparser? Or maybe just make this configurable from the + # command line with a global default? + # But the same font size might not be appropriate for all plots. + plt.rcParams["font.size"] = 20 + fig = plt.figure() + + conversion_list = list(sorted(set(df["X"]))) + pressure_list = list(sorted(set(df["P"]))) + n_conversion = len(conversion_list) + n_pressure = len(pressure_list) + error_lookup = { + (df["X"][i], df["P"][i]): df["surrogate-error"][i] + for i in range(len(df)) + } + error_array = [ + [error_lookup[conversion_list[j], pressure_list[i]] for i in range(n_pressure)] + for j in range(n_conversion) + ] + + cbar = legend + # "Maximum relative error" (Between reactor model and surrogate model, over + # all output values, for the intput values at which the optimization stopped) + cbar_kws = dict(label="Maximum relative error") + ax = sns.heatmap( + error_array, + cbar=cbar, + vmin=0.0, + vmax=0.05, + cbar_kws=cbar_kws, + square=True, + cmap="Reds", + ) + # Change "background color" used for missing data (where either full-space + # or surrogate model failed) + + if legend: + w, h = fig.get_size_inches() + # Getting the legend to appear in the right location is really annoying... + fig.set_size_inches(1.5*w, h) + legend_handles = [Patch(color="black", label="Failed simulation")] + ax.legend( + handles=legend_handles, + #ncol=1, + handlelength=0.75, + bbox_to_anchor=(2.5, 1.0), + #loc="upper right", + #borderaxespad=0, + ) + + xtick_positions = [i+0.5 for i in range(n_pressure)] + ytick_positions = [i+0.5 for i in range(n_conversion)] + xtick_labels = ["%1.2f" % (p / 1e6) if i%2 else "" for i, p in enumerate(pressure_list)] + ytick_labels = ["%1.2f" % x if i%2 else "" for i, x in enumerate(conversion_list)] + ax.set_xticks(xtick_positions, labels=xtick_labels) + ax.set_yticks(ytick_positions, labels=ytick_labels, rotation=0) + ax.set_xlabel("Pressure (MPa)") + ax.set_ylabel("Conversion") + + plt.gca().invert_yaxis() + ax.set_facecolor("black") + + return fig, ax + + +def main(args): + df = pd.read_csv(args.error_fpath) + + fig, ax = plot_surrogate_error( + df, + legend=not args.no_legend, + show_training_bounds=args.show_training_bounds, + ) + + if args.title is not None: + ax.set_title(args.title) + + fig.tight_layout() + + if not args.no_save: + if args.plot_fname is None: + plot_fname = os.path.basename(args.error_fpath) + data_ext = "." + plot_fname.split(".")[-1] + ext_len = len(data_ext) + plot_fname = plot_fname[:-ext_len] + ".pdf" + else: + plot_fname = args.plot_fname + + plot_fpath = os.path.join(args.results_dir, plot_fname) + fig.savefig(plot_fpath, transparent=not args.opaque) + + if args.show: + plt.show() + + +if __name__ == "__main__": + argparser = config.get_plot_argparser() + argparser.add_argument( + "error_fpath", help="Path to CSV file containing surrogate errors to plot" + ) + args = argparser.parse_args() + main(args) diff --git a/svi/auto_thermal_reformer/reactor_model.py b/svi/auto_thermal_reformer/reactor_model.py index 3b63d5f..ed51450 100644 --- a/svi/auto_thermal_reformer/reactor_model.py +++ b/svi/auto_thermal_reformer/reactor_model.py @@ -45,6 +45,7 @@ def add_reactor_model( flow_mol_h2o=300, flow_mol_gas=750, temp_gas=750, + initialize=True, ): """Add reactor unit model to provided Pyomo/IDAES model @@ -134,7 +135,8 @@ def add_reactor_model( m.fs.reformer_mix.initialize() propagate_state(arc=m.fs.connect) - m.fs.reformer.initialize() + if initialize: + m.fs.reformer.initialize() ######### SET REFORMER OUTLET PRESSURE ######### m.fs.reformer.outlet.pressure[0].fix(137895) @@ -158,6 +160,7 @@ def create_instance( flow_mol_h2o=300, flow_mol_gas=750, temp_gas=750, + initialize=True, ): # Set up global information m = ConcreteModel() @@ -172,6 +175,7 @@ def create_instance( flow_mol_h2o=flow_mol_h2o, flow_mol_gas=flow_mol_gas, temp_gas=temp_gas, + initialize=initialize, ) return m @@ -190,3 +194,29 @@ def create_instance( solver.solve(m, tee=True) tol = 1e-5 validate_solution(m, tolerance=tol) + + print("Variables") + print("---------") + for var in igraph.variables: + print(var.name) + print() + + print("Constraints") + print("-----------") + for con in igraph.constraints: + print(con.name) + print() + + vblocks, cblocks = igraph.block_triangularize() + for i, (vb, cb) in enumerate(zip(vblocks, cblocks)): + print(f"Block {i}") + print(f"=========") + print("Variables") + print("---------") + for var in vb: + print(f" {var.name}") + print("Constraints") + print("-----------") + for con in cb: + print(f" {con.name}") + print() diff --git a/svi/auto_thermal_reformer/results/parity_train_atr_NN.pdf b/svi/auto_thermal_reformer/results/parity_train_atr_NN.pdf deleted file mode 100644 index 4c16c14..0000000 Binary files a/svi/auto_thermal_reformer/results/parity_train_atr_NN.pdf and /dev/null differ diff --git a/svi/auto_thermal_reformer/results/parity_val_atr_nn.pdf b/svi/auto_thermal_reformer/results/parity_val_atr_nn.pdf deleted file mode 100644 index 5e04d86..0000000 Binary files a/svi/auto_thermal_reformer/results/parity_val_atr_nn.pdf and /dev/null differ diff --git a/svi/auto_thermal_reformer/run_alamo_sweep.py b/svi/auto_thermal_reformer/run_alamo_sweep.py index b83dca9..e4ab8f8 100644 --- a/svi/auto_thermal_reformer/run_alamo_sweep.py +++ b/svi/auto_thermal_reformer/run_alamo_sweep.py @@ -34,6 +34,7 @@ DEFAULT_SURROGATE_FNAME, ) import svi.auto_thermal_reformer.config as config +from svi.auto_thermal_reformer.compute_surrogate_error import compute_surrogate_error INVALID = None @@ -47,6 +48,11 @@ def main(): default="alamo-sweep.csv", help="Base file name for parameter sweep results", ) + argparser.add_argument( + "--compute-surrogate-error", + action="store_true", + help="Compute surrogate error (at termination point) and write to a separate results file", + ) args = argparser.parse_args() # TODO: This should be configurable by CLI @@ -55,6 +61,9 @@ def main(): df = {key: [] for key in config.PARAM_SWEEP_KEYS} + if args.compute_surrogate_error: + surrogate_error = {"X": [], "P": [], "surrogate-error": []} + """ The optimization problem to solve is the following: Maximize H2 composition in the product stream such that its minimum flow is 3500 mol/s, @@ -104,6 +113,28 @@ def main(): for key in config.PARAM_SWEEP_KEYS: if key not in ("X", "P", "Termination"): df[key].append(INVALID) + + # Regardless of whether the solve was successful, if we didn't encounter + # an error, attempt to compute the surrogate error at the point of + # termination. + if args.compute_surrogate_error: + # Computing this surrogate error involves two simulations. What + # if either fails? + surrogate_error["X"].append(X) + surrogate_error["P"].append(P) + try: + surr_err = compute_surrogate_error(m) + # Maximum relative error, defined as: + # |a - b| / max(|a|, |b|, 1) + max_surr_err = max(surr_err.values()) + surrogate_error["surrogate-error"].append(max_surr_err) + except ValueError: + # We failed to simulate either the full-space or surrogate + # reactor model. + surrogate_error["surrogate-error"].append(None) + + #finally: + # pass except ValueError: df[list(df.keys())[0]].append(X) df[list(df.keys())[1]].append(P) @@ -112,11 +143,34 @@ def main(): if key not in ("X", "P", "Termination"): df[key].append(INVALID) + if args.compute_surrogate_error: + # If we encountered an error, do not attempt to compute surrogate + # error. A solution has not been loaded, so the value we would + # get has nothing to do with the formulation or parameter values. + surrogate_error["X"].append(X) + surrogate_error["P"].append(P) + surrogate_error["surrogate-error"].append(INVALID) + df = pd.DataFrame(df) print(df) if not args.no_save: df.to_csv(output_fpath) + if args.compute_surrogate_error: + surrogate_error_df = pd.DataFrame(surrogate_error) + + # TODO: Allow specification of a file for surrogate error + sweep_basename = args.fname + if "." in sweep_basename: + experiment_extension = "." + sweep_basename.split(".")[-1] + name = sweep_basename[:-len(experiment_extension)] + error_fname = name + "-surrogate-error" + experiment_extension + else: + error_fname = sweep_basename + "-surrogate-error" + error_fpath = os.path.join(args.data_dir, error_fname) + + surrogate_error_df.to_csv(error_fpath) + if __name__ == "__main__": main() diff --git a/svi/auto_thermal_reformer/run_fullspace_sweep.py b/svi/auto_thermal_reformer/run_fullspace_sweep.py index 541a00b..2ba9f96 100644 --- a/svi/auto_thermal_reformer/run_fullspace_sweep.py +++ b/svi/auto_thermal_reformer/run_fullspace_sweep.py @@ -27,18 +27,21 @@ make_optimization_model, make_simulation_model, ) +from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP import svi.auto_thermal_reformer.config as config import pandas as pd import numpy as np - -df = {key: [] for key in config.PARAM_SWEEP_KEYS} - - INVALID = None +def calculate_condition_number(m): + nlp = PyomoNLP(m) + jac = nlp.evaluate_jacobian() + cond_num = np.linalg.cond(jac.toarray()) + return cond_num -def main(X,P): +def main(X,P, calc_condition_number = False): + m = make_optimization_model(X,P) # For instance 13 in the param sweep, these options give a quite interesting local @@ -50,6 +53,10 @@ def main(X,P): timer.tic("starting timer") print(f"Solving sample with X={X}, P={P}") results = solver.solve(m, tee=True, timer=htimer) + if calc_condition_number: + cond_num = calculate_condition_number(m) + else: + cond_num = None dT = timer.toc("end timer") f_eval_time = htimer.timers["solve"].timers["function"].total_time j_eval_time = htimer.timers["solve"].timers["jacobian"].total_time @@ -57,6 +64,8 @@ def main(X,P): df[list(df.keys())[0]].append(X) df[list(df.keys())[1]].append(P) df[list(df.keys())[2]].append(results.solver.termination_condition) + if new_key: + df[list(df.keys())[3]].append(cond_num) if pyo.check_optimal_termination(results): df["Time"].append(dT) df["Objective"].append(pyo.value(m.fs.product.mole_frac_comp[0,'H2'])) @@ -70,19 +79,37 @@ def main(X,P): else: # If the solver didn't converge, we don't care about the solve time, # the objective, or any of the degree of freedom values. - for key in config.PARAM_SWEEP_KEYS: - if key not in ("X", "P", "Termination"): + for key in DF_KEYS: + if key not in ("X", "P", "Termination", "Condition Number"): df[key].append(INVALID) if __name__ == "__main__": + argparser = config.get_sweep_argparser() + argparser.add_argument( "--fname", default="fullspace-sweep.csv", help="Base file name for parameter sweep results" ) + + argparser.add_argument( + "--calc-condition-number", + action="store_true", + help="Whether to calculate the condition number or not" + ) + args = argparser.parse_args() + + df = {key: [] for key in config.PARAM_SWEEP_KEYS} + index_to_insert = 3 + new_key = "Condition Number" if args.calc_condition_number else None + DF_KEYS = list(df.keys()) + if new_key: + DF_KEYS.insert(index_to_insert, new_key) + df = {key: [] for key in DF_KEYS} + xp_samples = config.get_parameter_samples(args) fpath = os.path.join(args.data_dir, args.fname) @@ -90,26 +117,26 @@ def main(X,P): for i, (X, P) in enumerate(xp_samples): print(f"Running sample {i} with X={X}, P={P}") try: - main(X,P) + main(X,P,args.calc_condition_number) except AssertionError: df[list(df.keys())[0]].append(X) df[list(df.keys())[1]].append(P) df[list(df.keys())[2]].append("AssertionError") - for key in config.PARAM_SWEEP_KEYS: + for key in DF_KEYS: if key not in ("X", "P", "Termination"): df[key].append(INVALID) except OverflowError: df[list(df.keys())[0]].append(X) df[list(df.keys())[1]].append(P) df[list(df.keys())[2]].append("OverflowError") - for key in config.PARAM_SWEEP_KEYS: + for key in DF_KEYS: if key not in ("X", "P", "Termination"): df[key].append(INVALID) except RuntimeError: df[list(df.keys())[0]].append(X) df[list(df.keys())[1]].append(P) df[list(df.keys())[2]].append("RuntimeError") - for key in config.PARAM_SWEEP_KEYS: + for key in DF_KEYS: if key not in ("X", "P", "Termination"): df[key].append(INVALID) diff --git a/svi/auto_thermal_reformer/run_nn_sweep.py b/svi/auto_thermal_reformer/run_nn_sweep.py index 3a006c1..723fb63 100644 --- a/svi/auto_thermal_reformer/run_nn_sweep.py +++ b/svi/auto_thermal_reformer/run_nn_sweep.py @@ -35,6 +35,7 @@ ) from idaes.core.surrogate.keras_surrogate import KerasSurrogate import svi.auto_thermal_reformer.config as config +from svi.auto_thermal_reformer.compute_surrogate_error import compute_surrogate_error INVALID = None @@ -70,6 +71,11 @@ def main(): " Must be 'full' or 'reduced'." ), ) + argparser.add_argument( + "--compute-surrogate-error", + action="store_true", + help="Compute surrogate error (at termination point) and write to a separate results file", + ) args = argparser.parse_args() @@ -77,12 +83,17 @@ def main(): if args.fname is None: sweep_fname = f"nn-sweep-{args.formulation}.csv" + else: + sweep_fname = args.fname surrogate_fname = os.path.join(args.data_dir, args.surrogate_fname) output_fpath = os.path.join(args.data_dir, sweep_fname) df = {key: [] for key in config.PARAM_SWEEP_KEYS} + if args.compute_surrogate_error: + surrogate_error = {"X": [], "P": [], "surrogate-error": []} + """ The optimization problem to solve is the following: Maximize H2 composition in the product stream such that its minimum flow is 3500 mol/s, @@ -134,6 +145,25 @@ def main(): for key in config.PARAM_SWEEP_KEYS: if key not in ("X", "P", "Termination"): df[key].append(INVALID) + + # Regardless of whether the solve was successful, if we didn't encounter + # an error, attempt to compute the surrogate error at the point of + # termination. + if args.compute_surrogate_error: + # Computing this surrogate error involves two simulations. What + # if either fails? + surrogate_error["X"].append(X) + surrogate_error["P"].append(P) + try: + surr_err = compute_surrogate_error(m) + # Maximum relative error, defined as: + # |a - b| / max(|a|, |b|, 1) + max_surr_err = max(surr_err.values()) + surrogate_error["surrogate-error"].append(max_surr_err) + except ValueError: + # We failed to simulate either the full-space or surrogate + # reactor model. + surrogate_error["surrogate-error"].append(None) except ValueError: df[list(df.keys())[0]].append(X) @@ -145,8 +175,31 @@ def main(): if key not in ("X", "P", "Termination"): df[key].append(INVALID) + if args.compute_surrogate_error: + # If we encountered an error, do not attempt to compute surrogate + # error. A solution has not been loaded, so the value we would + # get has nothing to do with the formulation or parameter values. + surrogate_error["X"].append(X) + surrogate_error["P"].append(P) + surrogate_error["surrogate-error"].append(INVALID) + df = pd.DataFrame(df) df.to_csv(output_fpath) + if args.compute_surrogate_error: + surrogate_error_df = pd.DataFrame(surrogate_error) + + # TODO: Allow specification of a file for surrogate error + sweep_basename = sweep_fname + if "." in sweep_basename: + experiment_extension = "." + sweep_basename.split(".")[-1] + name = sweep_basename[:-len(experiment_extension)] + error_fname = name + "-surrogate-error" + experiment_extension + else: + error_fname = sweep_basename + "-surrogate-error" + error_fpath = os.path.join(args.data_dir, error_fname) + + surrogate_error_df.to_csv(error_fpath) + if __name__ == "__main__": main() diff --git a/svi/cyipopt.py b/svi/cyipopt.py index be5f3a3..9504e44 100644 --- a/svi/cyipopt.py +++ b/svi/cyipopt.py @@ -14,6 +14,10 @@ from pyomo.core.base import Block, Objective, minimize from pyomo.opt import SolverStatus, SolverResults, TerminationCondition, ProblemSense from pyomo.opt.results.solution import Solution +import numpy as np +from scipy import sparse +from pyomo.contrib.incidence_analysis.triangularize import block_triangularize +from svi.nlp import get_gradient_of_lagrangian pyomo_nlp = attempt_import("pyomo.contrib.pynumero.interfaces.pyomo_nlp")[0] pyomo_grey_box = attempt_import("pyomo.contrib.pynumero.interfaces.pyomo_grey_box_nlp")[ @@ -38,6 +42,10 @@ def __init__(self, nlp, **kwds): timer = HierarchicalTimer() self._timer = timer + # We want this to be true so we can use 13argcallback below, but we can't + # use Pyomo v6.8.0 because of a version incompatibility with IDAES v2.4. + self._use_13arg_callback = True + super().__init__(nlp, **kwds) def solve(self, x, lagrange=None, zl=None, zu=None): @@ -297,9 +305,7 @@ def __init__(self): def __call__( self, nlp, - # Don't include this argument, for compatibility with current CyIpopt - # interface - #ipopt_problem, + ipopt_problem, alg_mod, iter_count, obj_value, @@ -327,3 +333,214 @@ def __call__( ls_trials, ) ) + +class ConditioningCallback: + + def __init__(self): + self.iterate_data = [] + self.condition_numbers = [] + + def __call__( + self, + nlp, + ipopt_problem, + alg_mod, + iter_count, + obj_value, + inf_pr, + inf_du, + mu, + d_norm, + regularization_size, + alpha_du, + alpha_pr, + ls_trials, + ): + self.iterate_data.append( + ( + alg_mod, + iter_count, + obj_value, + inf_pr, + inf_du, + mu, + d_norm, + regularization_size, + alpha_du, + alpha_pr, + ls_trials, + ) + ) + jac = nlp.evaluate_jacobian() + cond = np.linalg.cond(jac.toarray()) + self.condition_numbers.append(cond) + + +def extract_submatrix(coo, rows, cols): + new_rows = [] + new_cols = [] + new_data = [] + row_old2new = {r: i for i, r in enumerate(rows)} + col_old2new = {c: i for i, c in enumerate(cols)} + for r, c, d in zip(coo.row, coo.col, coo.data): + if r in row_old2new and c in col_old2new: + new_rows.append(row_old2new[r]) + new_cols.append(col_old2new[c]) + new_data.append(d) + new_coo = sparse.coo_matrix( + (new_data, (new_rows, new_cols)), + shape=(len(rows), len(cols)), + ) + return new_coo + + +class FullStateCallback: + + def __init__( + self, + include_condition=False, + include_block_condition=False, + dof_varnames=None, + ): + # TODO: We really don't want to re-use any any of this data from + # solve to solve. We need to store it, somehow, so we can access + # if after the callback function. Probably just need an explicit + # method to clear the cached data? Or check the NLP and make sure + # we're not being called with different NLPs? + # + # TODO: Option to compute condition numbers? + self.iterate_data = { + "alg_mod": [], + "iter_count": [], + "obj_value": [], + "inf_pr": [], + "inf_du": [], + "mu": [], + "d_norm": [], + "regularization_size": [], + "alpha_du": [], + "alpha_pr": [], + "ls_trials": [], + } + self.primal_values = {} + self.primal_residuals = {} + # TODO: Dual variables/residuals? + + self._cached_var_names = None + self._cached_con_names = None + + if include_condition: + self.condition_numbers = [] + else: + self.condition_numbers = None + + # Don't initialize with an empty list as we don't know how many elements + # this list will need. + self._include_block_condition = include_block_condition + if self._include_block_condition: + if dof_varnames is None: + raise RuntimeError("Degrees of freedom must be provided") + self._dof_varnames = dof_varnames + self.block_coords = None + self.block_condition_numbers = None + + def __call__( + self, + nlp, + ipopt_problem, + alg_mod, + iter_count, + obj_value, + inf_pr, + inf_du, + mu, + d_norm, + regularization_size, + alpha_du, + alpha_pr, + ls_trials, + ): + self.iterate_data["alg_mod"].append(alg_mod) + self.iterate_data["iter_count"].append(iter_count) + self.iterate_data["obj_value"].append(obj_value) + self.iterate_data["inf_pr"].append(inf_pr) + self.iterate_data["inf_du"].append(inf_du) + self.iterate_data["mu"].append(mu) + self.iterate_data["d_norm"].append(d_norm) + self.iterate_data["regularization_size"].append(regularization_size) + self.iterate_data["alpha_du"].append(alpha_du) + self.iterate_data["alpha_pr"].append(alpha_pr) + self.iterate_data["ls_trials"].append(ls_trials) + + iterate = ipopt_problem.get_current_iterate(scaled=False) + # NOTE: I'm assuming that the Ipopt iterate vector has the same order as the NLP. + # I'm pretty sure I've verified this before... + primal_values = iterate["x"] + if self._cached_var_names is None: + #self._cached_var_names = [var.name for var in nlp.get_pyomo_variables()] + self._cached_var_names = nlp.primals_names() + if not self.primal_values: + # If we haven't initialized this dict with variable names + for name, value in zip(self._cached_var_names, primal_values): + self.primal_values[name] = [value] + else: + # We have initialized with variable names + for name, value in zip(self._cached_var_names, primal_values): + self.primal_values[name].append(value) + + infeas = ipopt_problem.get_current_violations(scaled=False) + primal_residuals = infeas["g_violation"] + if self._cached_con_names is None: + #self._cached_con_names = [con.name for con in nlp.get_pyomo_constraints()] + self._cached_con_names = nlp.constraint_names() + if not self.primal_residuals: + for name, value in zip(self._cached_con_names, primal_residuals): + self.primal_residuals[name] = [value] + else: + for name, value in zip(self._cached_con_names, primal_residuals): + self.primal_residuals[name].append(value) + + # TODO: Could also track multipliers and dual infeasibility. This could be + # useful. + if self.condition_numbers is not None: + jac = nlp.evaluate_jacobian() + cond = np.linalg.cond(jac.toarray()) + self.condition_numbers.append(cond) + + if self._include_block_condition: + dof_varnames = set(self._dof_varnames) + var_coords = [i for i, vname in enumerate(self._cached_var_names) if vname not in dof_varnames] + square_varnames = [vname for i, vname in enumerate(self._cached_var_names) if vname not in dof_varnames] + name_to_new_coord = { + vname: i + for i, vname in enumerate(square_varnames) + if vname not in dof_varnames + } + con_coords = list(range(nlp.n_eq_constraints())) + vcoord_set = set(var_coords) + + jac = nlp.evaluate_jacobian_eq() + + row = [] + col = [] + data = [] + for r, c, val in zip(jac.row, jac.col, jac.data): + if c in vcoord_set: + new_coord = name_to_new_coord[self._cached_var_names[c]] + row.append(r) + col.append(new_coord) + data.append(val) + square_jac = sparse.coo_matrix((data, (row, col)), shape=(len(con_coords), len(var_coords))) + # Row and column blocks + rblocks, cblocks = block_triangularize(square_jac) + + if self.block_condition_numbers is None: + self.block_condition_numbers = { + f"block-{i}-cond": [] + for i in range(len(rblocks)) if len(rblocks[i]) > 1 + } + for i, (rb, cb) in enumerate(zip(rblocks, cblocks)): + if len(rb) > 1: + submat = extract_submatrix(square_jac, rb, cb) + cond = np.linalg.cond(submat.toarray()) + self.block_condition_numbers[f"block-{i}-cond"].append(cond) diff --git a/svi/nlp.py b/svi/nlp.py new file mode 100644 index 0000000..b459ae7 --- /dev/null +++ b/svi/nlp.py @@ -0,0 +1,95 @@ +from pyomo.common.collections import ComponentSet +import numpy as np +import scipy.sparse as sps + + +def get_reduced_hessian(nlp, dof_vars, lbmult, ubmult): + dofvarset = ComponentSet(dof_vars) + other_vars = [var for var in nlp.get_pyomo_variables() if var not in dofvarset] + eqcons = nlp.get_pyomo_equality_constraints() + dof_jac = nlp.extract_submatrix_jacobian(dof_vars, eqcons) + other_jac = nlp.extract_submatrix_jacobian(other_vars, eqcons) + ndof = len(dof_vars) + dof_nullbasis = np.identity(ndof) + lu = sps.linalg.splu(other_jac) + other_nullbasis = - lu.solve(dof_jac.toarray()) + nullbasis = np.vstack((dof_nullbasis, other_nullbasis)) + + varorder = dof_vars + other_vars + varindices = nlp.get_primal_indices(varorder) + + ineqcons = nlp.get_pyomo_inequality_constraints() + ineq_jac = nlp.extract_submatrix_jacobian(varorder, ineqcons) + ineq_val = np.diag(nlp.evaluate_ineq_constraints()) + ineq_val_inv = np.diag(1 / nlp.evaluate_ineq_constraints()) + ineq_duals = np.diag(nlp.get_duals_ineq()) + # This term is in the right order + ineq_term = ineq_jac.transpose() @ ineq_val_inv @ ineq_duals @ ineq_jac + + lbmult = np.diag(lbmult) + ubmult = np.diag(ubmult) + # TODO: Have to make sure values and multipliers are in the right order + primal_val = nlp.get_primals()[varindices] + primal_lb = nlp.primals_lb()[varindices] + primal_ub = nlp.primals_ub()[varindices] + # Have some variables equal to their (lower) bounds. These must not participate + # in the optimization. + active_lbs = np.where(primal_lb == primal_val)[0] + active_ubs = np.where(primal_ub == primal_val)[0] + primal_lb[active_lbs] -= 1e-8 + primal_ub[active_ubs] += 1e-8 + primal_ubslack_inv = np.diag(1 / (primal_ub - primal_val)) + primal_lbslack_inv = np.diag(1 / (primal_val - primal_lb)) + # TODO: Are these the right signs? + bound_term = lbmult @ primal_ubslack_inv + ubmult @ primal_lbslack_inv + + # Get Hessian in the proper order + # I'm just assuming that mutipliers here are using the right convention. + # This appears to be correct based on my grad-lag calculation. + hess = nlp.extract_submatrix_hessian_lag(varorder, varorder) + # Add terms for inequalities and bounds + hess += - ineq_term + bound_term + + rh = nullbasis.transpose() @ hess @ nullbasis + return rh + + +def project_onto(mat, coords): + # Schur complement of a dense matrix wrt the complement of coords + assert mat.shape[0] == mat.shape[1] + all_coords = np.arange(mat.shape[0]) + other_coords = all_coords[~np.isin(all_coords, coords)] + # Is this the right way to slice this matrix? + mat_11 = mat[coords, :][:, coords] + mat_22 = mat[other_coords, :][:, other_coords] + mat_12 = mat[coords, :][:, other_coords] + mat_21 = mat[other_coords, :][:, coords] + proj = mat_11 - mat_12.dot(np.linalg.solve(mat_22, mat_21)) + return proj + + +def get_gradient_of_lagrangian( + nlp, + primal_lb_multipliers, + primal_ub_multipliers, +): + # PyNumero NLPs contain constraint multipliers, but does not define a convention. + # We still need: + # - primal LB/UB multipliers + # We should not need slack multipliers (Ipopt should take care of this...) + grad_obj = nlp.evaluate_grad_objective() + + # There is no way this works. We will probably need to separate equality and + # inequality multipliers. + jac = nlp.evaluate_jacobian() + duals = nlp.get_duals() + # Each constraint gradient times its multiplier + conjac_term = jac.transpose().dot(duals) + + grad_lag = ( + - grad_obj + - conjac_term + + primal_lb_multipliers + - primal_ub_multipliers + ) + return grad_lag