diff --git a/doc/OnlineDocs/contributed_packages/pyros.rst b/doc/OnlineDocs/contributed_packages/pyros.rst index 4ef57fbf26c..23ec60f2e20 100644 --- a/doc/OnlineDocs/contributed_packages/pyros.rst +++ b/doc/OnlineDocs/contributed_packages/pyros.rst @@ -538,12 +538,16 @@ correspond to first-stage degrees of freedom. ... solve_master_globally=True, ... load_solution=False, ... ) - =========================================================================================== - PyROS: Pyomo Robust Optimization Solver ... - =========================================================================================== + ============================================================================== + PyROS: The Pyomo Robust Optimization Solver. ... - INFO: Robust optimal solution identified. Exiting PyROS. - + ------------------------------------------------------------------------------ + Robust optimal solution identified. + ------------------------------------------------------------------------------ + ... + ------------------------------------------------------------------------------ + All done. Exiting PyROS. + ============================================================================== >>> # === Query results === >>> time = results_1.time >>> iterations = results_1.iterations @@ -604,6 +608,8 @@ optional keyword argument ``decision_rule_order`` to the PyROS In this example, we select affine decision rules by setting ``decision_rule_order=1``: +.. _example-two-stg: + .. doctest:: :skipif: not (baron.available() and baron.license_is_valid()) @@ -625,11 +631,16 @@ In this example, we select affine decision rules by setting ... solve_master_globally=True, ... decision_rule_order=1, ... ) - =========================================================================================== - PyROS: Pyomo Robust Optimization Solver ... + ============================================================================== + PyROS: The Pyomo Robust Optimization Solver... ... - INFO: Robust optimal solution identified. Exiting PyROS. - + ------------------------------------------------------------------------------ + Robust optimal solution identified. + ------------------------------------------------------------------------------ + ... + ------------------------------------------------------------------------------ + All done. Exiting PyROS. + ============================================================================== >>> # === Compare final objective to the single-stage solution >>> two_stage_final_objective = round( ... pyo.value(results_2.final_objective_value), @@ -733,7 +744,289 @@ set size on the robust optimal objective function value and demonstrates the ease of implementing a price of robustness study for a given optimization problem under uncertainty. -.. note:: +PyROS Solver Log Output +------------------------------- + +The PyROS solver log output is controlled through the optional +``progress_logger`` argument, itself cast to +a standard Python logger (:py:class:`logging.Logger`) object +at the outset of a :meth:`~pyomo.contrib.pyros.PyROS.solve` call. +The level of detail of the solver log output +can be adjusted by adjusting the level of the +logger object; see :ref:`the following table `. +Note that by default, ``progress_logger`` is cast to a logger of level +:py:obj:`logging.INFO`. + +We refer the reader to the +:doc:`official Python logging library documentation ` +for customization of Python logger objects; +for a basic tutorial, see the :doc:`logging HOWTO `. + +.. _table-logging-levels: + +.. list-table:: PyROS solver log output at the various standard Python :py:mod:`logging` levels. + :widths: 10 50 + :header-rows: 1 + + * - Logging Level + - Output Messages + * - :py:obj:`logging.ERROR` + - * Information on the subproblem for which an exception was raised + by a subordinate solver + * Details about failure of the PyROS coefficient matching routine + * - :py:obj:`logging.WARNING` + - * Information about a subproblem not solved to an acceptable status + by the user-provided subordinate optimizers + * Invocation of a backup solver for a particular subproblem + * Caution about solution robustness guarantees in event that + user passes ``bypass_global_separation=True`` + * - :py:obj:`logging.INFO` + - * PyROS version, author, and disclaimer information + * Summary of user options + * Breakdown of model component statistics + * Iteration log table + * Termination details: message, timing breakdown, summary of statistics + * - :py:obj:`logging.DEBUG` + - * Termination outcomes and summary of statistics for + every master feasility, master, and DR polishing problem + * Progress updates for the separation procedure + * Separation subproblem initial point infeasibilities + * Summary of separation loop outcomes: performance constraints + violated, uncertain parameter scenario added to the + master problem + * Uncertain parameter scenarios added to the master problem + thus far + +An example of an output log produced through the default PyROS +progress logger is shown in +:ref:`the snippet that follows `. +Observe that the log contains the following information: + + +* **Introductory information** (lines 1--18). + Includes the version number, author + information, (UTC) time at which the solver was invoked, + and, if available, information on the local Git branch and + commit hash. +* **Summary of solver options** (lines 19--38). +* **Preprocessing information** (lines 39--41). + Wall time required for preprocessing + the deterministic model and associated components, + i.e. standardizing model components and adding the decision rule + variables and equations. +* **Model component statistics** (lines 42--58). + Breakdown of model component statistics. + Includes components added by PyROS, such as the decision rule variables + and equations. +* **Iteration log table** (lines 59--69). + Summary information on the problem iterates and subproblem outcomes. + The constituent columns are defined in detail in + :ref:`the table following the snippet `. +* **Termination message** (lines 70--71). Very brief summary of the termination outcome. +* **Timing statistics** (lines 72--88). + Tabulated breakdown of the solver timing statistics, based on a + :class:`pyomo.common.timing.HierarchicalTimer` printout. + The identifiers are as follows: + + * ``main``: Total time elapsed by the solver. + * ``main.dr_polishing``: Total time elapsed by the subordinate solvers + on polishing of the decision rules. + * ``main.global_separation``: Total time elapsed by the subordinate solvers + on global separation subproblems. + * ``main.local_separation``: Total time elapsed by the subordinate solvers + on local separation subproblems. + * ``main.master``: Total time elapsed by the subordinate solvers on + the master problems. + * ``main.master_feasibility``: Total time elapsed by the subordinate solvers + on the master feasibility problems. + * ``main.preprocessing``: Total preprocessing time. + * ``main.other``: Total overhead time. + +* **Termination statistics** (lines 89--94). Summary of statistics related to the + iterate at which PyROS terminates. +* **Exit message** (lines 95--96). + + +.. _solver-log-snippet: + +.. code-block:: text + :caption: PyROS solver output log for the :ref:`two-stage problem example `. + :linenos: + + ============================================================================== + PyROS: The Pyomo Robust Optimization Solver, v1.2.8. + Pyomo version: 6.7.0 + Commit hash: unknown + Invoked at UTC 2023-11-03T04:27:42.954101 + + Developed by: Natalie M. Isenberg (1), Jason A. F. Sherman (1), + John D. Siirola (2), Chrysanthos E. Gounaris (1) + (1) Carnegie Mellon University, Department of Chemical Engineering + (2) Sandia National Laboratories, Center for Computing Research + + The developers gratefully acknowledge support from the U.S. Department + of Energy's Institute for the Design of Advanced Energy Systems (IDAES). + ============================================================================== + ================================= DISCLAIMER ================================= + PyROS is still under development. + Please provide feedback and/or report any issues by creating a ticket at + https://github.com/Pyomo/pyomo/issues/new/choose + ============================================================================== + Solver options: + time_limit=None + keepfiles=False + tee=False + load_solution=True + objective_focus= + nominal_uncertain_param_vals=[0.13248000000000001, 4.97, 4.97, 1800] + decision_rule_order=1 + solve_master_globally=True + max_iter=-1 + robust_feasibility_tolerance=0.0001 + separation_priority_order={} + progress_logger= + backup_local_solvers=[] + backup_global_solvers=[] + subproblem_file_directory=None + bypass_local_separation=False + bypass_global_separation=False + p_robustness={} + ------------------------------------------------------------------------------ + Preprocessing... + Done preprocessing; required wall time of 0.175s. + ------------------------------------------------------------------------------ + Model statistics: + Number of variables : 62 + Epigraph variable : 1 + First-stage variables : 7 + Second-stage variables : 6 + State variables : 18 + Decision rule variables : 30 + Number of uncertain parameters : 4 + Number of constraints : 81 + Equality constraints : 24 + Coefficient matching constraints : 0 + Decision rule equations : 6 + All other equality constraints : 18 + Inequality constraints : 57 + First-stage inequalities (incl. certain var bounds) : 10 + Performance constraints (incl. var bounds) : 47 + ------------------------------------------------------------------------------ + Itn Objective 1-Stg Shift DR Shift #CViol Max Viol Wall Time (s) + ------------------------------------------------------------------------------ + 0 3.5838e+07 - - 5 1.8832e+04 1.198 + 1 3.5838e+07 7.4506e-09 1.6105e+03 7 3.7766e+04 2.893 + 2 3.6116e+07 2.7803e+05 1.2918e+03 8 1.3466e+06 4.732 + 3 3.6285e+07 1.6957e+05 5.8386e+03 6 4.8734e+03 6.740 + 4 3.6285e+07 1.4901e-08 3.3097e+03 1 3.5036e+01 9.099 + 5 3.6285e+07 2.9786e-10 3.3597e+03 6 2.9103e+00 11.588 + 6 3.6285e+07 7.4506e-07 8.7228e+02 5 4.1726e-01 14.360 + 7 3.6285e+07 7.4506e-07 8.1995e+02 0 9.3279e-10g 21.597 + ------------------------------------------------------------------------------ + Robust optimal solution identified. + ------------------------------------------------------------------------------ + Timing breakdown: + + Identifier ncalls cumtime percall % + ----------------------------------------------------------- + main 1 21.598 21.598 100.0 + ------------------------------------------------------ + dr_polishing 7 1.502 0.215 7.0 + global_separation 47 1.300 0.028 6.0 + local_separation 376 9.779 0.026 45.3 + master 8 5.385 0.673 24.9 + master_feasibility 7 0.531 0.076 2.5 + preprocessing 1 0.175 0.175 0.8 + other n/a 2.926 n/a 13.5 + ====================================================== + =========================================================== + + ------------------------------------------------------------------------------ + Termination stats: + Iterations : 8 + Solve time (wall s) : 21.598 + Final objective value : 3.6285e+07 + Termination condition : pyrosTerminationCondition.robust_optimal + ------------------------------------------------------------------------------ + All done. Exiting PyROS. + ============================================================================== + + +The iteration log table is designed to provide, in a concise manner, +important information about the progress of the iterative algorithm for +the problem of interest. +The constituent columns are defined in the +:ref:`table that follows `. + +.. _table-iteration-log-columns: + +.. list-table:: PyROS iteration log table columns. + :widths: 10 50 + :header-rows: 1 - Please provide feedback and/or report any problems by opening an issue on - the `Pyomo GitHub page `_. + * - Column Name + - Definition + * - Itn + - Iteration number. + * - Objective + - Master solution objective function value. + If the objective of the deterministic model provided + has a maximization sense, + then the negative of the objective function value is displayed. + Expect this value to trend upward as the iteration number + increases. + If the master problems are solved globally + (by passing ``solve_master_globally=True``), + then after the iteration number exceeds the number of uncertain parameters, + this value should be monotonically nondecreasing + as the iteration number is increased. + A dash ("-") is produced in lieu of a value if the master + problem of the current iteration is not solved successfully. + * - 1-Stg Shift + - Infinity norm of the difference between the first-stage + variable vectors of the master solutions of the current + and previous iterations. Expect this value to trend + downward as the iteration number increases. + A dash ("-") is produced in lieu of a value + if the current iteration number is 0, + there are no first-stage variables, + or the master problem of the current iteration is not solved successfully. + * - DR Shift + - Infinity norm of the difference between the decision rule + variable vectors of the master solutions of the current + and previous iterations. + Expect this value to trend downward as the iteration number increases. + An asterisk ("*") is appended to this value if the decision rules are + not successfully polished. + A dash ("-") is produced in lieu of a value + if the current iteration number is 0, + there are no decision rule variables, + or the master problem of the current iteration is not solved successfully. + * - #CViol + - Number of performance constraints found to be violated during + the separation step of the current iteration. + Unless a custom prioritization of the model's performance constraints + is specified (through the ``separation_priority_order`` argument), + expect this number to trend downward as the iteration number increases. + A "+" is appended if not all of the separation problems + were solved successfully, either due to custom prioritization, a time out, + or an issue encountered by the subordinate optimizers. + A dash ("-") is produced in lieu of a value if the separation + routine is not invoked during the current iteration. + * - Max Viol + - Maximum scaled performance constraint violation. + Expect this value to trend downward as the iteration number increases. + A 'g' is appended to the value if the separation problems were solved + globally during the current iteration. + A dash ("-") is produced in lieu of a value if the separation + routine is not invoked during the current iteration, or if there are + no performance constraints. + * - Wall time (s) + - Total time elapsed by the solver, in seconds, up to the end of the + current iteration. + + +Feedback and Reporting Issues +------------------------------- +Please provide feedback and/or report any problems by opening an issue on +the `Pyomo GitHub page `_. diff --git a/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py b/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py index cf1b8a2de23..e2d172f34f6 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py @@ -19,20 +19,20 @@ def main(): # Vars to estimate - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] # Data file_dirname = dirname(abspath(str(__file__))) - file_name = abspath(join(file_dirname, 'reactor_data.csv')) + file_name = abspath(join(file_dirname, "reactor_data.csv")) data = pd.read_csv(file_name) # Sum of squared error function def SSE(model, data): expr = ( - (float(data['ca']) - model.ca) ** 2 - + (float(data['cb']) - model.cb) ** 2 - + (float(data['cc']) - model.cc) ** 2 - + (float(data['cd']) - model.cd) ** 2 + (float(data.iloc[0]["ca"]) - model.ca) ** 2 + + (float(data.iloc[0]["cb"]) - model.cb) ** 2 + + (float(data.iloc[0]["cc"]) - model.cc) ** 2 + + (float(data.iloc[0]["cd"]) - model.cd) ** 2 ) return expr @@ -46,13 +46,13 @@ def SSE(model, data): bootstrap_theta = pest.theta_est_bootstrap(50) # Plot results - parmest.graphics.pairwise_plot(bootstrap_theta, title='Bootstrap theta') + parmest.graphics.pairwise_plot(bootstrap_theta, title="Bootstrap theta") parmest.graphics.pairwise_plot( bootstrap_theta, theta, 0.8, - ['MVN', 'KDE', 'Rect'], - title='Bootstrap theta with confidence regions', + ["MVN", "KDE", "Rect"], + title="Bootstrap theta with confidence regions", ) diff --git a/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py b/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py index 811571e20ed..cfd3891c00e 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py @@ -39,16 +39,16 @@ def generate_data(): data = pd.DataFrame() ndata = 200 # Normal distribution, mean = 3400, std = 500 - data['ca'] = 500 * np.random.randn(ndata) + 3400 + data["ca"] = 500 * np.random.randn(ndata) + 3400 # Random distribution between 500 and 1500 - data['cb'] = np.random.rand(ndata) * 1000 + 500 + data["cb"] = np.random.rand(ndata) * 1000 + 500 # Lognormal distribution - data['cc'] = np.random.lognormal(np.log(1600), 0.25, ndata) + data["cc"] = np.random.lognormal(np.log(1600), 0.25, ndata) # Triangular distribution between 1000 and 2000 - data['cd'] = np.random.triangular(1000, 1800, 3000, size=ndata) + data["cd"] = np.random.triangular(1000, 1800, 3000, size=ndata) - data['sv'] = sv_real - data['caf'] = caf_real + data["sv"] = sv_real + data["caf"] = caf_real return data @@ -61,10 +61,10 @@ def main(): # Define sum of squared error objective function for data rec def SSE(model, data): expr = ( - ((float(data['ca']) - model.ca) / float(data_std['ca'])) ** 2 - + ((float(data['cb']) - model.cb) / float(data_std['cb'])) ** 2 - + ((float(data['cc']) - model.cc) / float(data_std['cc'])) ** 2 - + ((float(data['cd']) - model.cd) / float(data_std['cd'])) ** 2 + ((float(data.iloc[0]["ca"]) - model.ca) / float(data_std["ca"])) ** 2 + + ((float(data.iloc[0]["cb"]) - model.cb) / float(data_std["cb"])) ** 2 + + ((float(data.iloc[0]["cc"]) - model.cc) / float(data_std["cc"])) ** 2 + + ((float(data.iloc[0]["cd"]) - model.cd) / float(data_std["cd"])) ** 2 ) return expr @@ -73,26 +73,26 @@ def SSE(model, data): pest = parmest.Estimator(reactor_design_model_for_datarec, data, theta_names, SSE) - obj, theta, data_rec = pest.theta_est(return_values=['ca', 'cb', 'cc', 'cd', 'caf']) + obj, theta, data_rec = pest.theta_est(return_values=["ca", "cb", "cc", "cd", "caf"]) print(obj) print(theta) parmest.graphics.grouped_boxplot( - data[['ca', 'cb', 'cc', 'cd']], - data_rec[['ca', 'cb', 'cc', 'cd']], - group_names=['Data', 'Data Rec'], + data[["ca", "cb", "cc", "cd"]], + data_rec[["ca", "cb", "cc", "cd"]], + group_names=["Data", "Data Rec"], ) ### Parameter estimation using reconciled data - theta_names = ['k1', 'k2', 'k3'] - data_rec['sv'] = data['sv'] + theta_names = ["k1", "k2", "k3"] + data_rec["sv"] = data["sv"] pest = parmest.Estimator(reactor_design_model, data_rec, theta_names, SSE) obj, theta = pest.theta_est() print(obj) print(theta) - theta_real = {'k1': 5.0 / 6.0, 'k2': 5.0 / 3.0, 'k3': 1.0 / 6000.0} + theta_real = {"k1": 5.0 / 6.0, "k2": 5.0 / 3.0, "k3": 1.0 / 6000.0} print(theta_real) diff --git a/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py b/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py index 95af53e63d3..6952a7fc733 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py @@ -20,11 +20,11 @@ def main(): # Vars to estimate - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] # Data file_dirname = dirname(abspath(str(__file__))) - file_name = abspath(join(file_dirname, 'reactor_data.csv')) + file_name = abspath(join(file_dirname, "reactor_data.csv")) data = pd.read_csv(file_name) # Create more data for the example @@ -37,10 +37,10 @@ def main(): # Sum of squared error function def SSE(model, data): expr = ( - (float(data['ca']) - model.ca) ** 2 - + (float(data['cb']) - model.cb) ** 2 - + (float(data['cc']) - model.cc) ** 2 - + (float(data['cd']) - model.cd) ** 2 + (float(data.iloc[0]["ca"]) - model.ca) ** 2 + + (float(data.iloc[0]["cb"]) - model.cb) ** 2 + + (float(data.iloc[0]["cc"]) - model.cc) ** 2 + + (float(data.iloc[0]["cd"]) - model.cd) ** 2 ) return expr @@ -68,7 +68,7 @@ def SSE(model, data): lNo = 25 lNo_samples = 5 bootstrap_samples = 20 - dist = 'MVN' + dist = "MVN" alphas = [0.7, 0.8, 0.9] results = pest.leaveNout_bootstrap_test( @@ -84,8 +84,8 @@ def SSE(model, data): bootstrap_results, theta_est_N, alpha, - ['MVN'], - title='Alpha: ' + str(alpha) + ', ' + str(theta_est_N.loc[0, alpha]), + ["MVN"], + title="Alpha: " + str(alpha) + ", " + str(theta_est_N.loc[0, alpha]), ) # Extract the percent of points that are within the alpha region diff --git a/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py b/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py index 13a40774740..a0fe6f22305 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py @@ -21,20 +21,20 @@ def main(): # Vars to estimate - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] # Data file_dirname = dirname(abspath(str(__file__))) - file_name = abspath(join(file_dirname, 'reactor_data.csv')) + file_name = abspath(join(file_dirname, "reactor_data.csv")) data = pd.read_csv(file_name) # Sum of squared error function def SSE(model, data): expr = ( - (float(data['ca']) - model.ca) ** 2 - + (float(data['cb']) - model.cb) ** 2 - + (float(data['cc']) - model.cc) ** 2 - + (float(data['cd']) - model.cd) ** 2 + (float(data.iloc[0]["ca"]) - model.ca) ** 2 + + (float(data.iloc[0]["cb"]) - model.cb) ** 2 + + (float(data.iloc[0]["cc"]) - model.cc) ** 2 + + (float(data.iloc[0]["cd"]) - model.cd) ** 2 ) return expr @@ -48,7 +48,7 @@ def SSE(model, data): k1 = [0.8, 0.85, 0.9] k2 = [1.6, 1.65, 1.7] k3 = [0.00016, 0.000165, 0.00017] - theta_vals = pd.DataFrame(list(product(k1, k2, k3)), columns=['k1', 'k2', 'k3']) + theta_vals = pd.DataFrame(list(product(k1, k2, k3)), columns=["k1", "k2", "k3"]) obj_at_theta = pest.objective_at_theta(theta_vals) # Run the likelihood ratio test @@ -56,7 +56,7 @@ def SSE(model, data): # Plot results parmest.graphics.pairwise_plot( - LR, theta, 0.9, title='LR results within 90% confidence region' + LR, theta, 0.9, title="LR results within 90% confidence region" ) diff --git a/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py b/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py index bc564cbdfd3..a92ac626fae 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py @@ -21,23 +21,23 @@ def main(): # Parameter estimation using multisensor data # Vars to estimate - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] # Data, includes multiple sensors for ca and cc file_dirname = dirname(abspath(str(__file__))) - file_name = abspath(join(file_dirname, 'reactor_data_multisensor.csv')) + file_name = abspath(join(file_dirname, "reactor_data_multisensor.csv")) data = pd.read_csv(file_name) # Sum of squared error function def SSE_multisensor(model, data): expr = ( - ((float(data['ca1']) - model.ca) ** 2) * (1 / 3) - + ((float(data['ca2']) - model.ca) ** 2) * (1 / 3) - + ((float(data['ca3']) - model.ca) ** 2) * (1 / 3) - + (float(data['cb']) - model.cb) ** 2 - + ((float(data['cc1']) - model.cc) ** 2) * (1 / 2) - + ((float(data['cc2']) - model.cc) ** 2) * (1 / 2) - + (float(data['cd']) - model.cd) ** 2 + ((float(data.iloc[0]["ca1"]) - model.ca) ** 2) * (1 / 3) + + ((float(data.iloc[0]["ca2"]) - model.ca) ** 2) * (1 / 3) + + ((float(data.iloc[0]["ca3"]) - model.ca) ** 2) * (1 / 3) + + (float(data.iloc[0]["cb"]) - model.cb) ** 2 + + ((float(data.iloc[0]["cc1"]) - model.cc) ** 2) * (1 / 2) + + ((float(data.iloc[0]["cc2"]) - model.cc) ** 2) * (1 / 2) + + (float(data.iloc[0]["cd"]) - model.cd) ** 2 ) return expr diff --git a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py b/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py index 334dfa264a4..581d3904c04 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py @@ -19,20 +19,20 @@ def main(): # Vars to estimate - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] # Data file_dirname = dirname(abspath(str(__file__))) - file_name = abspath(join(file_dirname, 'reactor_data.csv')) + file_name = abspath(join(file_dirname, "reactor_data.csv")) data = pd.read_csv(file_name) # Sum of squared error function def SSE(model, data): expr = ( - (float(data['ca']) - model.ca) ** 2 - + (float(data['cb']) - model.cb) ** 2 - + (float(data['cc']) - model.cc) ** 2 - + (float(data['cd']) - model.cd) ** 2 + (float(data.iloc[0]["ca"]) - model.ca) ** 2 + + (float(data.iloc[0]["cb"]) - model.cb) ** 2 + + (float(data.iloc[0]["cc"]) - model.cc) ** 2 + + (float(data.iloc[0]["cd"]) - model.cd) ** 2 ) return expr @@ -46,11 +46,11 @@ def SSE(model, data): k1_expected = 5.0 / 6.0 k2_expected = 5.0 / 3.0 k3_expected = 1.0 / 6000.0 - relative_error = abs(theta['k1'] - k1_expected) / k1_expected + relative_error = abs(theta["k1"] - k1_expected) / k1_expected assert relative_error < 0.05 - relative_error = abs(theta['k2'] - k2_expected) / k2_expected + relative_error = abs(theta["k2"] - k2_expected) / k2_expected assert relative_error < 0.05 - relative_error = abs(theta['k3'] - k3_expected) / k3_expected + relative_error = abs(theta["k3"] - k3_expected) / k3_expected assert relative_error < 0.05 diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py index f4cd6c8dbf5..16f65e236eb 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py +++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py @@ -37,10 +37,20 @@ def reactor_design_model(data): ) # m^3/(gmol min) # Inlet concentration of A, gmol/m^3 - model.caf = Param(initialize=float(data['caf']), within=PositiveReals) + if isinstance(data, dict) or isinstance(data, pd.Series): + model.caf = Param(initialize=float(data["caf"]), within=PositiveReals) + elif isinstance(data, pd.DataFrame): + model.caf = Param(initialize=float(data.iloc[0]["caf"]), within=PositiveReals) + else: + raise ValueError("Unrecognized data type.") # Space velocity (flowrate/volume) - model.sv = Param(initialize=float(data['sv']), within=PositiveReals) + if isinstance(data, dict) or isinstance(data, pd.Series): + model.sv = Param(initialize=float(data["sv"]), within=PositiveReals) + elif isinstance(data, pd.DataFrame): + model.sv = Param(initialize=float(data.iloc[0]["sv"]), within=PositiveReals) + else: + raise ValueError("Unrecognized data type.") # Outlet concentration of each component model.ca = Var(initialize=5000.0, within=PositiveReals) @@ -81,12 +91,12 @@ def main(): sv_values = [1.0 + v * 0.05 for v in range(1, 20)] caf = 10000 for sv in sv_values: - model = reactor_design_model({'caf': caf, 'sv': sv}) - solver = SolverFactory('ipopt') + model = reactor_design_model(pd.DataFrame(data={"caf": [caf], "sv": [sv]})) + solver = SolverFactory("ipopt") solver.solve(model) results.append([sv, caf, model.ca(), model.cb(), model.cc(), model.cd()]) - results = pd.DataFrame(results, columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd']) + results = pd.DataFrame(results, columns=["sv", "caf", "ca", "cb", "cc", "cd"]) print(results) diff --git a/pyomo/contrib/parmest/examples/semibatch/obj_at_theta.csv b/pyomo/contrib/parmest/examples/semibatch/obj_at_theta.csv new file mode 100644 index 00000000000..79f03e07dcd --- /dev/null +++ b/pyomo/contrib/parmest/examples/semibatch/obj_at_theta.csv @@ -0,0 +1,1009 @@ +,k1,k2,E1,E2,obj +0,4,40,29000,38000,667.4023645794207 +1,4,40,29000,38500,665.8312183437167 +2,4,40,29000,39000,672.7539769993407 +3,4,40,29000,39500,684.9503752463216 +4,4,40,29000,40000,699.985589093255 +5,4,40,29000,40500,716.1241770970677 +6,4,40,29000,41000,732.2023201586336 +7,4,40,29000,41500,747.4931745925483 +8,4,40,29500,38000,907.4405527163311 +9,4,40,29500,38500,904.2229271927299 +10,4,40,29500,39000,907.6942345285257 +11,4,40,29500,39500,915.4570013614677 +12,4,40,29500,40000,925.65401444575 +13,4,40,29500,40500,936.9348578520337 +14,4,40,29500,41000,948.3759339765711 +15,4,40,29500,41500,959.386491783636 +16,4,40,30000,38000,1169.8685711377334 +17,4,40,30000,38500,1166.2211505723928 +18,4,40,30000,39000,1167.702295374574 +19,4,40,30000,39500,1172.5517020611685 +20,4,40,30000,40000,1179.3820406408263 +21,4,40,30000,40500,1187.1698633839655 +22,4,40,30000,41000,1195.2047840919602 +23,4,40,30000,41500,1203.0241101248102 +24,4,40,30500,38000,1445.9591944684807 +25,4,40,30500,38500,1442.6632745483 +26,4,40,30500,39000,1443.1982444457385 +27,4,40,30500,39500,1446.2833842279929 +28,4,40,30500,40000,1450.9012120934779 +29,4,40,30500,40500,1456.295140290636 +30,4,40,30500,41000,1461.9350767569827 +31,4,40,30500,41500,1467.4715014446226 +32,4,40,31000,38000,1726.8744994061449 +33,4,40,31000,38500,1724.2679845375048 +34,4,40,31000,39000,1724.4550886870552 +35,4,40,31000,39500,1726.5124587129135 +36,4,40,31000,40000,1729.7061680616455 +37,4,40,31000,40500,1733.48893482641 +38,4,40,31000,41000,1737.4753558920438 +39,4,40,31000,41500,1741.4093763605517 +40,4,40,31500,38000,2004.1978135112938 +41,4,40,31500,38500,2002.2807839860222 +42,4,40,31500,39000,2002.3676405166086 +43,4,40,31500,39500,2003.797808439923 +44,4,40,31500,40000,2006.048051591001 +45,4,40,31500,40500,2008.7281679153625 +46,4,40,31500,41000,2011.5626384878237 +47,4,40,31500,41500,2014.3675286347284 +48,4,80,29000,38000,845.8197358579285 +49,4,80,29000,38500,763.5039795545781 +50,4,80,29000,39000,709.8529964173656 +51,4,80,29000,39500,679.4215539491266 +52,4,80,29000,40000,666.4876088521157 +53,4,80,29000,40500,665.978271760966 +54,4,80,29000,41000,673.7240200504901 +55,4,80,29000,41500,686.4763909417914 +56,4,80,29500,38000,1042.519415429413 +57,4,80,29500,38500,982.8097210678039 +58,4,80,29500,39000,942.2990207573541 +59,4,80,29500,39500,917.9550916645245 +60,4,80,29500,40000,906.3116029967189 +61,4,80,29500,40500,904.0326666308792 +62,4,80,29500,41000,908.1964630052729 +63,4,80,29500,41500,916.4222043837499 +64,4,80,30000,38000,1271.1030403496538 +65,4,80,30000,38500,1227.7527550544085 +66,4,80,30000,39000,1197.433957624904 +67,4,80,30000,39500,1178.447676126182 +68,4,80,30000,40000,1168.645219243497 +69,4,80,30000,40500,1165.7995210546096 +70,4,80,30000,41000,1167.8586496250396 +71,4,80,30000,41500,1173.0949214020527 +72,4,80,30500,38000,1520.8220402652044 +73,4,80,30500,38500,1489.2563260709424 +74,4,80,30500,39000,1466.8099189128857 +75,4,80,30500,39500,1452.4352624958806 +76,4,80,30500,40000,1444.7074679423818 +77,4,80,30500,40500,1442.0820578624343 +78,4,80,30500,41000,1443.099006489627 +79,4,80,30500,41500,1446.5106517200784 +80,4,80,31000,38000,1781.149136032395 +81,4,80,31000,38500,1758.2414369536502 +82,4,80,31000,39000,1741.891639711003 +83,4,80,31000,39500,1731.358661496594 +84,4,80,31000,40000,1725.6231647999593 +85,4,80,31000,40500,1723.5757174297378 +86,4,80,31000,41000,1724.1680229486278 +87,4,80,31000,41500,1726.5050840601884 +88,4,80,31500,38000,2042.8335948845602 +89,4,80,31500,38500,2026.3067503042414 +90,4,80,31500,39000,2014.5720701940838 +91,4,80,31500,39500,2007.0463766643977 +92,4,80,31500,40000,2002.9647983728314 +93,4,80,31500,40500,2001.5163951989875 +94,4,80,31500,41000,2001.9474217001339 +95,4,80,31500,41500,2003.6204088755821 +96,4,120,29000,38000,1176.0713512305115 +97,4,120,29000,38500,1016.8213383282462 +98,4,120,29000,39000,886.0136231565133 +99,4,120,29000,39500,789.0101180066036 +100,4,120,29000,40000,724.5420056133441 +101,4,120,29000,40500,686.6877602625062 +102,4,120,29000,41000,668.8129085873959 +103,4,120,29000,41500,665.1167761036883 +104,4,120,29500,38000,1263.887274509128 +105,4,120,29500,38500,1155.6528408872423 +106,4,120,29500,39000,1066.393539894248 +107,4,120,29500,39500,998.9931006471243 +108,4,120,29500,40000,952.36314487701 +109,4,120,29500,40500,923.4000293372077 +110,4,120,29500,41000,908.407361383214 +111,4,120,29500,41500,903.8136176328255 +112,4,120,30000,38000,1421.1418235449091 +113,4,120,30000,38500,1347.114022652679 +114,4,120,30000,39000,1285.686103704643 +115,4,120,30000,39500,1238.2456448658272 +116,4,120,30000,40000,1204.3526810790904 +117,4,120,30000,40500,1182.4272879027071 +118,4,120,30000,41000,1170.3447810121902 +119,4,120,30000,41500,1165.8422968073423 +120,4,120,30500,38000,1625.5588911535713 +121,4,120,30500,38500,1573.5546642859429 +122,4,120,30500,39000,1530.1592840718379 +123,4,120,30500,39500,1496.2087139473604 +124,4,120,30500,40000,1471.525855239756 +125,4,120,30500,40500,1455.2084749904016 +126,4,120,30500,41000,1445.9160840082027 +127,4,120,30500,41500,1442.1255377330835 +128,4,120,31000,38000,1855.8467211183756 +129,4,120,31000,38500,1818.4368412235558 +130,4,120,31000,39000,1787.25956706785 +131,4,120,31000,39500,1762.8169908546402 +132,4,120,31000,40000,1744.9825741661596 +133,4,120,31000,40500,1733.136625016882 +134,4,120,31000,41000,1726.3352245899828 +135,4,120,31000,41500,1723.492199933745 +136,4,120,31500,38000,2096.6479813687533 +137,4,120,31500,38500,2069.3606691038876 +138,4,120,31500,39000,2046.792043575205 +139,4,120,31500,39500,2029.2128703900223 +140,4,120,31500,40000,2016.4664599897606 +141,4,120,31500,40500,2008.054814885348 +142,4,120,31500,41000,2003.2622557140814 +143,4,120,31500,41500,2001.289784483679 +144,7,40,29000,38000,149.32898706737052 +145,7,40,29000,38500,161.04814413969586 +146,7,40,29000,39000,187.87801343005242 +147,7,40,29000,39500,223.00789161520424 +148,7,40,29000,40000,261.66779887964003 +149,7,40,29000,40500,300.676316191238 +150,7,40,29000,41000,338.04021206995765 +151,7,40,29000,41500,372.6191631389286 +152,7,40,29500,38000,276.6495061185777 +153,7,40,29500,38500,282.1304583501965 +154,7,40,29500,39000,300.91417483065254 +155,7,40,29500,39500,327.24304394350395 +156,7,40,29500,40000,357.0561976596432 +157,7,40,29500,40500,387.61662064170207 +158,7,40,29500,41000,417.1836349752378 +159,7,40,29500,41500,444.73705844573243 +160,7,40,30000,38000,448.0380830353589 +161,7,40,30000,38500,448.8094536459122 +162,7,40,30000,39000,460.77530593327293 +163,7,40,30000,39500,479.342874472736 +164,7,40,30000,40000,501.20694459059405 +165,7,40,30000,40500,524.0971649678811 +166,7,40,30000,41000,546.539334134893 +167,7,40,30000,41500,567.6447156158981 +168,7,40,30500,38000,657.9909416906933 +169,7,40,30500,38500,655.7465129488842 +170,7,40,30500,39000,662.5420970804985 +171,7,40,30500,39500,674.8914651553109 +172,7,40,30500,40000,690.2111920703564 +173,7,40,30500,40500,706.6833639709198 +174,7,40,30500,41000,723.0994507096715 +175,7,40,30500,41500,738.7096013891406 +176,7,40,31000,38000,899.1769906655776 +177,7,40,31000,38500,895.4391505892945 +178,7,40,31000,39000,898.7695629120826 +179,7,40,31000,39500,906.603316771593 +180,7,40,31000,40000,916.9811481373996 +181,7,40,31000,40500,928.4913367709245 +182,7,40,31000,41000,940.1744934710283 +183,7,40,31000,41500,951.4199286075984 +184,7,40,31500,38000,1163.093373675207 +185,7,40,31500,38500,1159.0457727559028 +186,7,40,31500,39000,1160.3831770028223 +187,7,40,31500,39500,1165.2451698296604 +188,7,40,31500,40000,1172.1768190340001 +189,7,40,31500,40500,1180.1105659428963 +190,7,40,31500,41000,1188.3083929833688 +191,7,40,31500,41500,1196.29112579565 +192,7,80,29000,38000,514.0332369183081 +193,7,80,29000,38500,329.3645784712966 +194,7,80,29000,39000,215.73000998706416 +195,7,80,29000,39500,162.37338399591852 +196,7,80,29000,40000,149.8401793263549 +197,7,80,29000,40500,162.96125998112578 +198,7,80,29000,41000,191.173279165834 +199,7,80,29000,41500,227.2781971491003 +200,7,80,29500,38000,623.559246695578 +201,7,80,29500,38500,448.60620511421484 +202,7,80,29500,39000,344.21940687907573 +203,7,80,29500,39500,292.9758707105001 +204,7,80,29500,40000,277.07670134364804 +205,7,80,29500,40500,283.5158840045542 +206,7,80,29500,41000,303.33951582820265 +207,7,80,29500,41500,330.43357046741954 +208,7,80,30000,38000,732.5907387079073 +209,7,80,30000,38500,593.1926567994672 +210,7,80,30000,39000,508.5638538704666 +211,7,80,30000,39500,464.47881763522037 +212,7,80,30000,40000,448.0394620671692 +213,7,80,30000,40500,449.64309860415494 +214,7,80,30000,41000,462.4490598612332 +215,7,80,30000,41500,481.6323506247537 +216,7,80,30500,38000,871.1163930229344 +217,7,80,30500,38500,771.1320563649375 +218,7,80,30500,39000,707.8872660015606 +219,7,80,30500,39500,672.6612145133173 +220,7,80,30500,40000,657.4974157809264 +221,7,80,30500,40500,656.0835852491216 +222,7,80,30500,41000,663.6006958125331 +223,7,80,30500,41500,676.460675405631 +224,7,80,31000,38000,1053.1852617390061 +225,7,80,31000,38500,984.3647109805877 +226,7,80,31000,39000,938.6158531749268 +227,7,80,31000,39500,911.4268280093535 +228,7,80,31000,40000,898.333365348419 +229,7,80,31000,40500,895.3996527486954 +230,7,80,31000,41000,899.3556288533885 +231,7,80,31000,41500,907.6180684887955 +232,7,80,31500,38000,1274.2255948763498 +233,7,80,31500,38500,1226.5236809533717 +234,7,80,31500,39000,1193.4538731398666 +235,7,80,31500,39500,1172.8105398345213 +236,7,80,31500,40000,1162.0692230240734 +237,7,80,31500,40500,1158.7461521476607 +238,7,80,31500,41000,1160.6173577210805 +239,7,80,31500,41500,1165.840315694716 +240,7,120,29000,38000,1325.2409732290193 +241,7,120,29000,38500,900.8063148840154 +242,7,120,29000,39000,629.9300352098937 +243,7,120,29000,39500,413.81648033893424 +244,7,120,29000,40000,257.3116751690404 +245,7,120,29000,40500,177.89217179438947 +246,7,120,29000,41000,151.58366848473491 +247,7,120,29000,41500,157.56967437251706 +248,7,120,29500,38000,1211.2807882170853 +249,7,120,29500,38500,956.936161969002 +250,7,120,29500,39000,753.3050086992201 +251,7,120,29500,39500,528.2452647799327 +252,7,120,29500,40000,382.62610532894917 +253,7,120,29500,40500,308.44199089882375 +254,7,120,29500,41000,280.3893024671524 +255,7,120,29500,41500,280.4028092582749 +256,7,120,30000,38000,1266.5740351143413 +257,7,120,30000,38500,1084.3028700477778 +258,7,120,30000,39000,834.2392498526193 +259,7,120,30000,39500,650.7560171314304 +260,7,120,30000,40000,537.7846910878052 +261,7,120,30000,40500,477.3001078155485 +262,7,120,30000,41000,451.6865380286754 +263,7,120,30000,41500,448.14911508024613 +264,7,120,30500,38000,1319.6603196780936 +265,7,120,30500,38500,1102.3027489012372 +266,7,120,30500,39000,931.2523583659847 +267,7,120,30500,39500,807.0833484596384 +268,7,120,30500,40000,727.4852710400268 +269,7,120,30500,40500,682.1437030344305 +270,7,120,30500,41000,660.7859329989657 +271,7,120,30500,41500,655.6001132492668 +272,7,120,31000,38000,1330.5306924865326 +273,7,120,31000,38500,1195.9190861202942 +274,7,120,31000,39000,1086.0328080422887 +275,7,120,31000,39500,1005.4160637517409 +276,7,120,31000,40000,951.2021706290612 +277,7,120,31000,40500,918.1457644271304 +278,7,120,31000,41000,901.0511005554887 +279,7,120,31000,41500,895.4599964465793 +280,7,120,31500,38000,1447.8365822059013 +281,7,120,31500,38500,1362.3417347939844 +282,7,120,31500,39000,1292.382727215108 +283,7,120,31500,39500,1239.1826828976662 +284,7,120,31500,40000,1201.6474412465277 +285,7,120,31500,40500,1177.5235955796813 +286,7,120,31500,41000,1164.1761722345295 +287,7,120,31500,41500,1158.9997785002718 +288,10,40,29000,38000,33.437068437082054 +289,10,40,29000,38500,58.471249815534996 +290,10,40,29000,39000,101.41937628542912 +291,10,40,29000,39500,153.80690200519626 +292,10,40,29000,40000,209.66451461551316 +293,10,40,29000,40500,265.03070792175197 +294,10,40,29000,41000,317.46079310177566 +295,10,40,29000,41500,365.59950388342645 +296,10,40,29500,38000,70.26818405688635 +297,10,40,29500,38500,87.96463718548947 +298,10,40,29500,39000,122.58188233160993 +299,10,40,29500,39500,166.2478945807132 +300,10,40,29500,40000,213.48669617414316 +301,10,40,29500,40500,260.67953961944477 +302,10,40,29500,41000,305.5877041218316 +303,10,40,29500,41500,346.95612213021155 +304,10,40,30000,38000,153.67588703371362 +305,10,40,30000,38500,164.07504103479005 +306,10,40,30000,39000,190.0800160661499 +307,10,40,30000,39500,224.61382980242837 +308,10,40,30000,40000,262.79232847382445 +309,10,40,30000,40500,301.38687703450415 +310,10,40,30000,41000,338.38536686093164 +311,10,40,30000,41500,372.6399011703545 +312,10,40,30500,38000,284.2936286531718 +313,10,40,30500,38500,288.4690608277705 +314,10,40,30500,39000,306.44667517621144 +315,10,40,30500,39500,332.20122250191986 +316,10,40,30500,40000,361.5566690083291 +317,10,40,30500,40500,391.72755224929614 +318,10,40,30500,41000,420.95317535960476 +319,10,40,30500,41500,448.2049230608669 +320,10,40,31000,38000,459.03140021766137 +321,10,40,31000,38500,458.71477027519967 +322,10,40,31000,39000,469.9910751800656 +323,10,40,31000,39500,488.05850105225426 +324,10,40,31000,40000,509.5204701455629 +325,10,40,31000,40500,532.0674969691778 +326,10,40,31000,41000,554.2088430693509 +327,10,40,31000,41500,575.0485839499048 +328,10,40,31500,38000,672.2476845983564 +329,10,40,31500,38500,669.2240508488649 +330,10,40,31500,39000,675.4956226836405 +331,10,40,31500,39500,687.447764319295 +332,10,40,31500,40000,702.4395430742891 +333,10,40,31500,40500,718.6279487347668 +334,10,40,31500,41000,734.793684592168 +335,10,40,31500,41500,750.1821072409286 +336,10,80,29000,38000,387.7617282731497 +337,10,80,29000,38500,195.33642612593002 +338,10,80,29000,39000,82.7306931465102 +339,10,80,29000,39500,35.13436471793541 +340,10,80,29000,40000,33.521138659248706 +341,10,80,29000,40500,61.47395975053128 +342,10,80,29000,41000,106.71403229340167 +343,10,80,29000,41500,160.56068704487473 +344,10,80,29500,38000,459.63404601804103 +345,10,80,29500,38500,258.7453720995899 +346,10,80,29500,39000,135.96435731320256 +347,10,80,29500,39500,80.2685095017944 +348,10,80,29500,40000,70.86302366453106 +349,10,80,29500,40500,90.43203026480438 +350,10,80,29500,41000,126.7844695901737 +351,10,80,29500,41500,171.63682876805044 +352,10,80,30000,38000,564.1463320344325 +353,10,80,30000,38500,360.75718124523866 +354,10,80,30000,39000,231.70119191254307 +355,10,80,30000,39500,170.74752201483128 +356,10,80,30000,40000,154.7149036950422 +357,10,80,30000,40500,166.10596450541493 +358,10,80,30000,41000,193.3351721194443 +359,10,80,30000,41500,228.78394172417038 +360,10,80,30500,38000,689.6797223218513 +361,10,80,30500,38500,484.8023695265838 +362,10,80,30500,39000,363.5979340028588 +363,10,80,30500,39500,304.67857102688225 +364,10,80,30500,40000,285.29210000833734 +365,10,80,30500,40500,290.0135917456113 +366,10,80,30500,41000,308.8672169492536 +367,10,80,30500,41500,335.3210332569182 +368,10,80,31000,38000,789.946106942773 +369,10,80,31000,38500,625.7722360026959 +370,10,80,31000,39000,528.6063264942235 +371,10,80,31000,39500,478.6863763478618 +372,10,80,31000,40000,459.5026243189753 +373,10,80,31000,40500,459.6982093164963 +374,10,80,31000,41000,471.6790024321937 +375,10,80,31000,41500,490.3034492109124 +376,10,80,31500,38000,912.3540488244158 +377,10,80,31500,38500,798.2135101409633 +378,10,80,31500,39000,727.746684419146 +379,10,80,31500,39500,689.0119464356724 +380,10,80,31500,40000,672.0757202772029 +381,10,80,31500,40500,669.678339553036 +382,10,80,31500,41000,676.5761221409929 +383,10,80,31500,41500,688.9934449650118 +384,10,120,29000,38000,1155.1165164624408 +385,10,120,29000,38500,840.2641727088946 +386,10,120,29000,39000,506.9102636732852 +387,10,120,29000,39500,265.5278912452038 +388,10,120,29000,40000,116.39516513179322 +389,10,120,29000,40500,45.2088092745619 +390,10,120,29000,41000,30.22267557153353 +391,10,120,29000,41500,51.06063746392809 +392,10,120,29500,38000,1343.7868459826054 +393,10,120,29500,38500,977.9852373227346 +394,10,120,29500,39000,594.632756549817 +395,10,120,29500,39500,346.2478773329187 +396,10,120,29500,40000,180.23082247413407 +397,10,120,29500,40500,95.81649989178923 +398,10,120,29500,41000,71.0837801649128 +399,10,120,29500,41500,82.84289818279714 +400,10,120,30000,38000,1532.9333545384934 +401,10,120,30000,38500,1012.2223350568845 +402,10,120,30000,39000,688.4884716222766 +403,10,120,30000,39500,464.6206903113392 +404,10,120,30000,40000,283.5644748300334 +405,10,120,30000,40500,190.27593217865416 +406,10,120,30000,41000,158.0192279691727 +407,10,120,30000,41500,161.3611926772337 +408,10,120,30500,38000,1349.3785399811063 +409,10,120,30500,38500,1014.785480110738 +410,10,120,30500,39000,843.0316833766408 +411,10,120,30500,39500,589.4543896730125 +412,10,120,30500,40000,412.3358512291996 +413,10,120,30500,40500,324.11715620464133 +414,10,120,30500,41000,290.17588242984766 +415,10,120,30500,41500,287.56857384673356 +416,10,120,31000,38000,1328.0973931040146 +417,10,120,31000,38500,1216.5659656437845 +418,10,120,31000,39000,928.4831767181619 +419,10,120,31000,39500,700.3115484040329 +420,10,120,31000,40000,565.0876352458171 +421,10,120,31000,40500,494.44016026435037 +422,10,120,31000,41000,464.38005437182983 +423,10,120,31000,41500,458.7614573733091 +424,10,120,31500,38000,1473.1154650008834 +425,10,120,31500,38500,1195.943614951571 +426,10,120,31500,39000,990.2486604382486 +427,10,120,31500,39500,843.1390407497395 +428,10,120,31500,40000,751.2746391170706 +429,10,120,31500,40500,700.215375503209 +430,10,120,31500,41000,676.1585052687219 +431,10,120,31500,41500,669.5907920932743 +432,13,40,29000,38000,49.96352152045025 +433,13,40,29000,38500,83.75104994958261 +434,13,40,29000,39000,136.8176091795391 +435,13,40,29000,39500,199.91486685466407 +436,13,40,29000,40000,266.4367154860076 +437,13,40,29000,40500,331.97224579940524 +438,13,40,29000,41000,393.8001583706036 +439,13,40,29000,41500,450.42425363084493 +440,13,40,29500,38000,29.775721038786923 +441,13,40,29500,38500,57.37673742631121 +442,13,40,29500,39000,103.49161398239501 +443,13,40,29500,39500,159.3058253852367 +444,13,40,29500,40000,218.60083223764073 +445,13,40,29500,40500,277.2507278183831 +446,13,40,29500,41000,332.7141278886951 +447,13,40,29500,41500,383.58832292300576 +448,13,40,30000,38000,47.72263852005472 +449,13,40,30000,38500,68.07581028940402 +450,13,40,30000,39000,106.13974628945516 +451,13,40,30000,39500,153.58449949683063 +452,13,40,30000,40000,204.62393623358633 +453,13,40,30000,40500,255.44513025602419 +454,13,40,30000,41000,303.69954914051766 +455,13,40,30000,41500,348.0803709720354 +456,13,40,30500,38000,110.9331168284094 +457,13,40,30500,38500,123.63361262704746 +458,13,40,30500,39000,153.02654433825705 +459,13,40,30500,39500,191.40769947472756 +460,13,40,30500,40000,233.503841403055 +461,13,40,30500,40500,275.8557790922913 +462,13,40,30500,41000,316.32529882763697 +463,13,40,30500,41500,353.7060432094809 +464,13,40,31000,38000,221.90608823073939 +465,13,40,31000,38500,227.67026441593657 +466,13,40,31000,39000,248.62107049869064 +467,13,40,31000,39500,277.9507605389158 +468,13,40,31000,40000,311.0267471957685 +469,13,40,31000,40500,344.8024031161673 +470,13,40,31000,41000,377.3761144228052 +471,13,40,31000,41500,407.6529635071056 +472,13,40,31500,38000,378.8738382757093 +473,13,40,31500,38500,379.39748335944216 +474,13,40,31500,39000,393.01223361732553 +475,13,40,31500,39500,414.10238059122855 +476,13,40,31500,40000,438.8024282436204 +477,13,40,31500,40500,464.5348067190265 +478,13,40,31500,41000,489.6621039898805 +479,13,40,31500,41500,513.2163939332803 +480,13,80,29000,38000,364.387588581215 +481,13,80,29000,38500,184.2902007673634 +482,13,80,29000,39000,81.57192155036655 +483,13,80,29000,39500,42.54811210095659 +484,13,80,29000,40000,49.897338772663076 +485,13,80,29000,40500,87.84229516509882 +486,13,80,29000,41000,143.85451969447664 +487,13,80,29000,41500,208.71467984917848 +488,13,80,29500,38000,382.5794635435733 +489,13,80,29500,38500,188.38619353711718 +490,13,80,29500,39000,75.75749359688277 +491,13,80,29500,39500,29.27891251986562 +492,13,80,29500,40000,29.794874961934568 +493,13,80,29500,40500,60.654888662698205 +494,13,80,29500,41000,109.25801388824325 +495,13,80,29500,41500,166.6311093454692 +496,13,80,30000,38000,448.97795526074816 +497,13,80,30000,38500,238.44530107604737 +498,13,80,30000,39000,112.34545890264337 +499,13,80,30000,39500,56.125871791222835 +500,13,80,30000,40000,48.29987461781518 +501,13,80,30000,40500,70.7900626637678 +502,13,80,30000,41000,110.76865376691964 +503,13,80,30000,41500,159.50197316936024 +504,13,80,30500,38000,547.7818730461195 +505,13,80,30500,38500,332.92604070423494 +506,13,80,30500,39000,193.80760050280742 +507,13,80,30500,39500,128.3457644087917 +508,13,80,30500,40000,112.23915895822442 +509,13,80,30500,40500,125.96369396512564 +510,13,80,30500,41000,156.67918617660013 +511,13,80,30500,41500,196.05195109523765 +512,13,80,31000,38000,682.8591931963246 +513,13,80,31000,38500,457.56562267948556 +514,13,80,31000,39000,313.6380169123524 +515,13,80,31000,39500,245.13531819580908 +516,13,80,31000,40000,223.54473391202873 +517,13,80,31000,40500,229.60752111202834 +518,13,80,31000,41000,251.42377424735136 +519,13,80,31000,41500,281.48720903016886 +520,13,80,31500,38000,807.925638050234 +521,13,80,31500,38500,588.686585641994 +522,13,80,31500,39000,464.0488586698228 +523,13,80,31500,39500,402.69214492641095 +524,13,80,31500,40000,380.13626165363934 +525,13,80,31500,40500,380.8064948609387 +526,13,80,31500,41000,395.05186915919086 +527,13,80,31500,41500,416.70193045600774 +528,13,120,29000,38000,1068.8279454397398 +529,13,120,29000,38500,743.0012805963486 +530,13,120,29000,39000,451.2538301167544 +531,13,120,29000,39500,235.4154251166075 +532,13,120,29000,40000,104.73720814447498 +533,13,120,29000,40500,46.91983990671749 +534,13,120,29000,41000,42.81092192562316 +535,13,120,29000,41500,74.33530639171506 +536,13,120,29500,38000,1133.1178848710972 +537,13,120,29500,38500,824.0745323788527 +538,13,120,29500,39000,499.10867111401996 +539,13,120,29500,39500,256.1626809904186 +540,13,120,29500,40000,107.68599585294751 +541,13,120,29500,40500,38.18533662516749 +542,13,120,29500,41000,25.499608203619154 +543,13,120,29500,41500,49.283537699300375 +544,13,120,30000,38000,1292.409871290162 +545,13,120,30000,38500,994.669572829704 +546,13,120,30000,39000,598.9783697712826 +547,13,120,30000,39500,327.47348408537925 +548,13,120,30000,40000,156.82634841081907 +549,13,120,30000,40500,71.30833688875883 +550,13,120,30000,41000,47.72389750130817 +551,13,120,30000,41500,62.1982461882982 +552,13,120,30500,38000,1585.8797221278146 +553,13,120,30500,38500,1144.66688416451 +554,13,120,30500,39000,692.6651441690645 +555,13,120,30500,39500,441.98837639874046 +556,13,120,30500,40000,251.56311435857728 +557,13,120,30500,40500,149.79670413140468 +558,13,120,30500,41000,115.52645596043719 +559,13,120,30500,41500,120.44019473389324 +560,13,120,31000,38000,1702.7625866892163 +561,13,120,31000,38500,1071.7854750250656 +562,13,120,31000,39000,807.8943299034604 +563,13,120,31000,39500,588.672223513561 +564,13,120,31000,40000,376.44658358671404 +565,13,120,31000,40500,269.2159719426485 +566,13,120,31000,41000,229.41660529009877 +567,13,120,31000,41500,226.78274707181976 +568,13,120,31500,38000,1331.3523701291767 +569,13,120,31500,38500,1151.2055268669133 +570,13,120,31500,39000,1006.811285091974 +571,13,120,31500,39500,702.0053094629535 +572,13,120,31500,40000,515.9081891614829 +573,13,120,31500,40500,423.8652275555525 +574,13,120,31500,41000,386.4939696097151 +575,13,120,31500,41500,379.8118453367429 +576,16,40,29000,38000,106.1025746852808 +577,16,40,29000,38500,145.32590128581407 +578,16,40,29000,39000,204.74804378224422 +579,16,40,29000,39500,274.6339266648551 +580,16,40,29000,40000,347.9667393938497 +581,16,40,29000,40500,420.03753452490974 +582,16,40,29000,41000,487.9353932879741 +583,16,40,29000,41500,550.0623063219693 +584,16,40,29500,38000,54.65040870471303 +585,16,40,29500,38500,88.94089091627293 +586,16,40,29500,39000,142.72223808288405 +587,16,40,29500,39500,206.63598763907422 +588,16,40,29500,40000,273.99851593521134 +589,16,40,29500,40500,340.34861536649436 +590,16,40,29500,41000,402.935270882596 +591,16,40,29500,41500,460.2471155081633 +592,16,40,30000,38000,29.788548081995298 +593,16,40,30000,38500,57.96323252610644 +594,16,40,30000,39000,104.92815906834525 +595,16,40,30000,39500,161.71867032726158 +596,16,40,30000,40000,222.01677586338877 +597,16,40,30000,40500,281.6349465235367 +598,16,40,30000,41000,337.99683241119567 +599,16,40,30000,41500,389.68271710858414 +600,16,40,30500,38000,42.06569536892785 +601,16,40,30500,38500,62.95145274276575 +602,16,40,30500,39000,101.93860830594608 +603,16,40,30500,39500,150.47910837525734 +604,16,40,30500,40000,202.65388851823258 +605,16,40,30500,40500,254.5724108541227 +606,16,40,30500,41000,303.84403622726694 +607,16,40,30500,41500,349.1422884543064 +608,16,40,31000,38000,99.21707896667829 +609,16,40,31000,38500,112.24153596941301 +610,16,40,31000,39000,142.5186177618655 +611,16,40,31000,39500,182.02836955332134 +612,16,40,31000,40000,225.3201896575212 +613,16,40,31000,40500,268.83705389232614 +614,16,40,31000,41000,310.3895932135811 +615,16,40,31000,41500,348.7480165565453 +616,16,40,31500,38000,204.30418825821732 +617,16,40,31500,38500,210.0759235359138 +618,16,40,31500,39000,231.7643258544752 +619,16,40,31500,39500,262.1512494310348 +620,16,40,31500,40000,296.3864127264238 +621,16,40,31500,40500,331.30743171999035 +622,16,40,31500,41000,364.95322314895554 +623,16,40,31500,41500,396.20142191205844 +624,16,80,29000,38000,399.5975649320935 +625,16,80,29000,38500,225.6318269911425 +626,16,80,29000,39000,127.97354075513151 +627,16,80,29000,39500,93.73584101549991 +628,16,80,29000,40000,106.43084032022394 +629,16,80,29000,40500,150.51245762256931 +630,16,80,29000,41000,213.24213500046466 +631,16,80,29000,41500,285.0426423013882 +632,16,80,29500,38000,371.37706087096393 +633,16,80,29500,38500,189.77150413822454 +634,16,80,29500,39000,86.22375488959844 +635,16,80,29500,39500,46.98714814001572 +636,16,80,29500,40000,54.596900621760675 +637,16,80,29500,40500,93.12033833747024 +638,16,80,29500,41000,149.89341227947025 +639,16,80,29500,41500,215.5937000584367 +640,16,80,30000,38000,388.43657991253195 +641,16,80,30000,38500,190.77121362008674 +642,16,80,30000,39000,76.28535232335287 +643,16,80,30000,39500,29.152860363695716 +644,16,80,30000,40000,29.820972887404942 +645,16,80,30000,40500,61.320203047752464 +646,16,80,30000,41000,110.82086782062603 +647,16,80,30000,41500,169.197767615573 +648,16,80,30500,38000,458.8964339917103 +649,16,80,30500,38500,239.547928886725 +650,16,80,30500,39000,109.02338779317503 +651,16,80,30500,39500,50.888746196140914 +652,16,80,30500,40000,42.73606982375976 +653,16,80,30500,40500,65.75935122724029 +654,16,80,30500,41000,106.68884313872147 +655,16,80,30500,41500,156.54100549486617 +656,16,80,31000,38000,561.7385153195615 +657,16,80,31000,38500,335.5692026144635 +658,16,80,31000,39000,188.0383015831574 +659,16,80,31000,39500,118.2318539104416 +660,16,80,31000,40000,100.81000168801492 +661,16,80,31000,40500,114.72014539486217 +662,16,80,31000,41000,146.2992492326178 +663,16,80,31000,41500,186.8074429488408 +664,16,80,31500,38000,697.9937997454152 +665,16,80,31500,38500,466.42234442578484 +666,16,80,31500,39000,306.52125608515166 +667,16,80,31500,39500,230.54692639209762 +668,16,80,31500,40000,206.461121102699 +669,16,80,31500,40500,212.23429887269359 +670,16,80,31500,41000,234.70913795495554 +671,16,80,31500,41500,265.8143069252357 +672,16,120,29000,38000,1085.688903883652 +673,16,120,29000,38500,750.2887000017752 +674,16,120,29000,39000,469.92662852990964 +675,16,120,29000,39500,267.1560282754928 +676,16,120,29000,40000,146.06299930062625 +677,16,120,29000,40500,95.28836772053619 +678,16,120,29000,41000,97.41466545178946 +679,16,120,29000,41500,135.3804131941845 +680,16,120,29500,38000,1079.5576154477903 +681,16,120,29500,38500,751.2932384998761 +682,16,120,29500,39000,458.27083477307207 +683,16,120,29500,39500,240.9658024131812 +684,16,120,29500,40000,109.3801465044384 +685,16,120,29500,40500,51.274139057659724 +686,16,120,29500,41000,47.36446629605638 +687,16,120,29500,41500,79.42944320845996 +688,16,120,30000,38000,1139.3792936518537 +689,16,120,30000,38500,833.7979589668842 +690,16,120,30000,39000,507.805443202025 +691,16,120,30000,39500,259.93892964607977 +692,16,120,30000,40000,108.7341499557062 +693,16,120,30000,40500,38.152937143498605 +694,16,120,30000,41000,25.403985123518716 +695,16,120,30000,41500,49.72822589160786 +696,16,120,30500,38000,1285.0396277304772 +697,16,120,30500,38500,1025.254169031627 +698,16,120,30500,39000,622.5890550779666 +699,16,120,30500,39500,333.3353043756717 +700,16,120,30500,40000,155.70268128051293 +701,16,120,30500,40500,66.84125446522368 +702,16,120,30500,41000,42.25187049753978 +703,16,120,30500,41500,56.98314898830595 +704,16,120,31000,38000,1595.7993459811262 +705,16,120,31000,38500,1252.8886556470425 +706,16,120,31000,39000,731.4408383874198 +707,16,120,31000,39500,451.0090473423308 +708,16,120,31000,40000,251.5086563526081 +709,16,120,31000,40500,141.8915050063955 +710,16,120,31000,41000,104.67474675582574 +711,16,120,31000,41500,109.1609567535697 +712,16,120,31500,38000,1942.3896021770768 +713,16,120,31500,38500,1197.207050908449 +714,16,120,31500,39000,812.6818768064074 +715,16,120,31500,39500,611.45532452889 +716,16,120,31500,40000,380.63642711770643 +717,16,120,31500,40500,258.5514125337487 +718,16,120,31500,41000,213.48518421250665 +719,16,120,31500,41500,209.58134396574906 +720,19,40,29000,38000,169.3907733115706 +721,19,40,29000,38500,212.23331960093145 +722,19,40,29000,39000,275.9376503672959 +723,19,40,29000,39500,350.4301397081139 +724,19,40,29000,40000,428.40863665493924 +725,19,40,29000,40500,504.955113902399 +726,19,40,29000,41000,577.023450987656 +727,19,40,29000,41500,642.9410032211753 +728,19,40,29500,38000,102.40889356493292 +729,19,40,29500,38500,141.19036226103668 +730,19,40,29500,39000,200.19333708701748 +731,19,40,29500,39500,269.6750686488757 +732,19,40,29500,40000,342.6217886299377 +733,19,40,29500,40500,414.33044375626207 +734,19,40,29500,41000,481.89521316730713 +735,19,40,29500,41500,543.7211700546151 +736,19,40,30000,38000,51.95330426445395 +737,19,40,30000,38500,85.69656829127965 +738,19,40,30000,39000,138.98376466247876 +739,19,40,30000,39500,202.43251598105033 +740,19,40,30000,40000,269.3557903452929 +741,19,40,30000,40500,335.2960133312316 +742,19,40,30000,41000,397.50658847538665 +743,19,40,30000,41500,454.47903112410967 +744,19,40,30500,38000,28.864802790801026 +745,19,40,30500,38500,56.32899754732796 +746,19,40,30500,39000,102.69825523352162 +747,19,40,30500,39500,158.95118263535466 +748,19,40,30500,40000,218.75241957992617 +749,19,40,30500,40500,277.9122290233915 +750,19,40,30500,41000,333.8561815041273 +751,19,40,30500,41500,385.1662652901447 +752,19,40,31000,38000,43.72359701781447 +753,19,40,31000,38500,63.683967347844224 +754,19,40,31000,39000,101.95579433282329 +755,19,40,31000,39500,149.8826019475827 +756,19,40,31000,40000,201.50605279789198 +757,19,40,31000,40500,252.92391570754876 +758,19,40,31000,41000,301.7431453727685 +759,19,40,31000,41500,346.6368192781496 +760,19,40,31500,38000,104.05710998615942 +761,19,40,31500,38500,115.95783594434451 +762,19,40,31500,39000,145.42181873662554 +763,19,40,31500,39500,184.26373455825217 +764,19,40,31500,40000,226.97066340897095 +765,19,40,31500,40500,269.96403356902357 +766,19,40,31500,41000,311.04753558871505 +767,19,40,31500,41500,348.98866332680115 +768,19,80,29000,38000,453.1314944429312 +769,19,80,29000,38500,281.24067760117225 +770,19,80,29000,39000,185.83730378881882 +771,19,80,29000,39500,154.25726305915472 +772,19,80,29000,40000,170.2912737797755 +773,19,80,29000,40500,218.38979299191152 +774,19,80,29000,41000,285.604024444273 +775,19,80,29000,41500,362.0858325427657 +776,19,80,29500,38000,400.06299682217264 +777,19,80,29500,38500,224.41725666435008 +778,19,80,29500,39000,125.58476107530382 +779,19,80,29500,39500,90.55733834394478 +780,19,80,29500,40000,102.67519971027264 +781,19,80,29500,40500,146.27807815967392 +782,19,80,29500,41000,208.57372904155937 +783,19,80,29500,41500,279.9669583078214 +784,19,80,30000,38000,376.1594584816549 +785,19,80,30000,38500,191.30452808298463 +786,19,80,30000,39000,85.63116084217559 +787,19,80,30000,39500,45.10487847849711 +788,19,80,30000,40000,51.88389644342952 +789,19,80,30000,40500,89.78942817703852 +790,19,80,30000,41000,146.0393555385696 +791,19,80,30000,41500,211.26567367707352 +792,19,80,30500,38000,401.874315275947 +793,19,80,30500,38500,197.55305366608133 +794,19,80,30500,39000,79.00348967857379 +795,19,80,30500,39500,29.602719961568614 +796,19,80,30500,40000,28.980451378502487 +797,19,80,30500,40500,59.63541802023186 +798,19,80,30500,41000,108.48607655362268 +799,19,80,30500,41500,166.30589286399507 +800,19,80,31000,38000,484.930958445979 +801,19,80,31000,38500,254.27552635537404 +802,19,80,31000,39000,116.75543721560439 +803,19,80,31000,39500,54.77547840250418 +804,19,80,31000,40000,44.637472658824976 +805,19,80,31000,40500,66.50466903927668 +806,19,80,31000,41000,106.62737262508298 +807,19,80,31000,41500,155.8310688191254 +808,19,80,31500,38000,595.6094306603337 +809,19,80,31500,38500,359.60040819463063 +810,19,80,31500,39000,201.85328967228585 +811,19,80,31500,39500,126.24442464793601 +812,19,80,31500,40000,106.07388975142673 +813,19,80,31500,40500,118.52358345403363 +814,19,80,31500,41000,149.1597537162607 +815,19,80,31500,41500,188.94964975523197 +816,19,120,29000,38000,1133.9213841599772 +817,19,120,29000,38500,793.9759807804692 +818,19,120,29000,39000,516.5580425563733 +819,19,120,29000,39500,318.60172051726147 +820,19,120,29000,40000,201.662212274693 +821,19,120,29000,40500,154.47522945829064 +822,19,120,29000,41000,160.28049502033574 +823,19,120,29000,41500,202.35345983501588 +824,19,120,29500,38000,1091.6343400395158 +825,19,120,29500,38500,754.9332443184217 +826,19,120,29500,39000,472.1777992591152 +827,19,120,29500,39500,267.03951846894995 +828,19,120,29500,40000,144.25558152688114 +829,19,120,29500,40500,92.40384156679512 +830,19,120,29500,41000,93.81833253459942 +831,19,120,29500,41500,131.24753560710644 +832,19,120,30000,38000,1092.719296892266 +833,19,120,30000,38500,764.7065490850255 +834,19,120,30000,39000,467.2268758064373 +835,19,120,30000,39500,244.9367732985332 +836,19,120,30000,40000,110.00996333393202 +837,19,120,30000,40500,49.96381544207811 +838,19,120,30000,41000,44.9298739569088 +839,19,120,30000,41500,76.25447129089613 +840,19,120,30500,38000,1160.6160120981158 +841,19,120,30500,38500,865.5953188304933 +842,19,120,30500,39000,531.1657093741892 +843,19,120,30500,39500,271.98520008106277 +844,19,120,30500,40000,114.03616090967407 +845,19,120,30500,40500,39.74252227099571 +846,19,120,30500,41000,25.07176465285551 +847,19,120,30500,41500,48.298794094852724 +848,19,120,31000,38000,1304.8870694342509 +849,19,120,31000,38500,1089.6854636757826 +850,19,120,31000,39000,668.6632735260521 +851,19,120,31000,39500,356.7751012890747 +852,19,120,31000,40000,168.32491564142487 +853,19,120,31000,40500,72.82648063377391 +854,19,120,31000,41000,45.02326687759286 +855,19,120,31000,41500,58.13111530831655 +856,19,120,31500,38000,1645.2697164013964 +857,19,120,31500,38500,1373.859712069864 +858,19,120,31500,39000,787.3948673670299 +859,19,120,31500,39500,483.60546305948367 +860,19,120,31500,40000,273.4285373433001 +861,19,120,31500,40500,153.21079535396908 +862,19,120,31500,41000,111.21299419905313 +863,19,120,31500,41500,113.52006337929113 +864,22,40,29000,38000,229.2032513971666 +865,22,40,29000,38500,274.65023153674116 +866,22,40,29000,39000,341.4424739822062 +867,22,40,29000,39500,419.2624324130753 +868,22,40,29000,40000,500.6022690006133 +869,22,40,29000,40500,580.3923016374031 +870,22,40,29000,41000,655.4874207991389 +871,22,40,29000,41500,724.1595537770351 +872,22,40,29500,38000,155.45206306046595 +873,22,40,29500,38500,197.41588482427002 +874,22,40,29500,39000,260.1641484982308 +875,22,40,29500,39500,333.666918810689 +876,22,40,29500,40000,410.66541588422854 +877,22,40,29500,40500,486.276072112155 +878,22,40,29500,41000,557.4760464927683 +879,22,40,29500,41500,622.6057687448293 +880,22,40,30000,38000,90.70026588811803 +881,22,40,30000,38500,128.41239603755494 +882,22,40,30000,39000,186.27261386900233 +883,22,40,30000,39500,254.5802373859711 +884,22,40,30000,40000,326.3686182341553 +885,22,40,30000,40500,396.9735001502319 +886,22,40,30000,41000,463.5155278718613 +887,22,40,30000,41500,524.414569320113 +888,22,40,30500,38000,44.551475763397946 +889,22,40,30500,38500,76.95264448905411 +890,22,40,30500,39000,128.85898727872572 +891,22,40,30500,39500,190.91422001003792 +892,22,40,30500,40000,256.4755613806196 +893,22,40,30500,40500,321.125224208803 +894,22,40,30500,41000,382.14434919800453 +895,22,40,30500,41500,438.03974322333033 +896,22,40,31000,38000,28.101321546315717 +897,22,40,31000,38500,53.867829756398805 +898,22,40,31000,39000,98.57619184859544 +899,22,40,31000,39500,153.19473192134507 +900,22,40,31000,40000,211.4202434313414 +901,22,40,31000,40500,269.09905982026265 +902,22,40,31000,41000,323.68306330754416 +903,22,40,31000,41500,373.76836451736045 +904,22,40,31500,38000,51.648288279447364 +905,22,40,31500,38500,69.56074881661863 +906,22,40,31500,39000,105.91402675097291 +907,22,40,31500,39500,151.99456204656389 +908,22,40,31500,40000,201.85995274525234 +909,22,40,31500,40500,251.63807959916412 +910,22,40,31500,41000,298.9593498669657 +911,22,40,31500,41500,342.50888994628025 +912,22,80,29000,38000,507.5440336860194 +913,22,80,29000,38500,336.42019672232965 +914,22,80,29000,39000,242.21016116765423 +915,22,80,29000,39500,212.33396533224905 +916,22,80,29000,40000,230.67632355958136 +917,22,80,29000,40500,281.6224662955561 +918,22,80,29000,41000,352.0457411487133 +919,22,80,29000,41500,431.89288175778637 +920,22,80,29500,38000,443.2889283037078 +921,22,80,29500,38500,270.0648237630224 +922,22,80,29500,39000,173.57666711629645 +923,22,80,29500,39500,141.06258420240613 +924,22,80,29500,40000,156.18412870159142 +925,22,80,29500,40500,203.33105261575707 +926,22,80,29500,41000,269.5552387411201 +927,22,80,29500,41500,345.03801326123767 +928,22,80,30000,38000,395.34177505602497 +929,22,80,30000,38500,217.11094192826982 +930,22,80,30000,39000,116.38535634181476 +931,22,80,30000,39500,79.94742924888467 +932,22,80,30000,40000,90.84706550421288 +933,22,80,30000,40500,133.26308067939766 +934,22,80,30000,41000,194.36064414396228 +935,22,80,30000,41500,264.56059537656466 +936,22,80,30500,38000,382.0341866812038 +937,22,80,30500,38500,191.65621311671836 +938,22,80,30500,39000,82.3318677587146 +939,22,80,30500,39500,39.44606931321677 +940,22,80,30500,40000,44.476166488763134 +941,22,80,30500,40500,80.84561981845566 +942,22,80,30500,41000,135.62459431793735 +943,22,80,30500,41500,199.42208168600175 +944,22,80,31000,38000,425.5181957619983 +945,22,80,31000,38500,210.2667219741389 +946,22,80,31000,39000,84.97041062888985 +947,22,80,31000,39500,31.593073529038755 +948,22,80,31000,40000,28.407154164211214 +949,22,80,31000,40500,57.05446633976857 +950,22,80,31000,41000,104.10423883907688 +951,22,80,31000,41500,160.23135976433713 +952,22,80,31500,38000,527.5015417150911 +953,22,80,31500,38500,282.29650611769665 +954,22,80,31500,39000,134.62881845323489 +955,22,80,31500,39500,66.62736532046851 +956,22,80,31500,40000,52.9918858786988 +957,22,80,31500,40500,72.36913743145999 +958,22,80,31500,41000,110.38003828747726 +959,22,80,31500,41500,157.65470091455973 +960,22,120,29000,38000,1186.823326813257 +961,22,120,29000,38500,844.3317816964005 +962,22,120,29000,39000,567.7367986440256 +963,22,120,29000,39500,371.79782508970567 +964,22,120,29000,40000,256.9261857702517 +965,22,120,29000,40500,211.85466060592006 +966,22,120,29000,41000,220.09534855737033 +967,22,120,29000,41500,265.02731793490034 +968,22,120,29500,38000,1128.4568915685559 +969,22,120,29500,38500,787.7709648712951 +970,22,120,29500,39000,508.4832626962424 +971,22,120,29500,39500,308.52654841064975 +972,22,120,29500,40000,190.01030358402707 +973,22,120,29500,40500,141.62663282114926 +974,22,120,29500,41000,146.40704203984612 +975,22,120,29500,41500,187.48734389188584 +976,22,120,30000,38000,1094.7007205604846 +977,22,120,30000,38500,757.7313528729464 +978,22,120,30000,39000,471.282561364766 +979,22,120,30000,39500,262.0412520036699 +980,22,120,30000,40000,136.26956239282435 +981,22,120,30000,40500,82.4268827471484 +982,22,120,30000,41000,82.3695177584498 +983,22,120,30000,41500,118.51210034475737 +984,22,120,30500,38000,1111.0872182758205 +985,22,120,30500,38500,787.2204655558988 +986,22,120,30500,39000,481.85960605002055 +987,22,120,30500,39500,250.28740868446397 +988,22,120,30500,40000,109.21968920710272 +989,22,120,30500,40500,45.51600269221681 +990,22,120,30500,41000,38.172157811051115 +991,22,120,30500,41500,67.73748641348168 +992,22,120,31000,38000,1193.3958874354898 +993,22,120,31000,38500,923.0731791194576 +994,22,120,31000,39000,573.4457650536078 +995,22,120,31000,39500,294.2980811757103 +996,22,120,31000,40000,124.86249624679849 +997,22,120,31000,40500,43.948524347749846 +998,22,120,31000,41000,25.582084045731808 +999,22,120,31000,41500,46.36268252714472 +1000,22,120,31500,38000,1336.0993444856913 +1001,22,120,31500,38500,1194.893001664831 +1002,22,120,31500,39000,740.6584250286721 +1003,22,120,31500,39500,397.18127104230757 +1004,22,120,31500,40000,194.20390582893873 +1005,22,120,31500,40500,88.22588964369922 +1006,22,120,31500,41000,54.97797247760634 +1007,22,120,31500,41500,64.88195101638016 diff --git a/pyomo/contrib/parmest/examples/semibatch/semibatch.py b/pyomo/contrib/parmest/examples/semibatch/semibatch.py index 8cda262c019..6762531a338 100644 --- a/pyomo/contrib/parmest/examples/semibatch/semibatch.py +++ b/pyomo/contrib/parmest/examples/semibatch/semibatch.py @@ -34,11 +34,20 @@ def generate_model(data): + # if data is a file name, then load file first + if isinstance(data, str): + file_name = data + try: + with open(file_name, "r") as infile: + data = json.load(infile) + except: + raise RuntimeError(f"Could not read {file_name} as json") + # unpack and fix the data - cameastemp = data['Ca_meas'] - cbmeastemp = data['Cb_meas'] - ccmeastemp = data['Cc_meas'] - trmeastemp = data['Tr_meas'] + cameastemp = data["Ca_meas"] + cbmeastemp = data["Cb_meas"] + ccmeastemp = data["Cc_meas"] + trmeastemp = data["Tr_meas"] cameas = {} cbmeas = {} @@ -79,9 +88,9 @@ def generate_model(data): m.Vc = Param(initialize=0.07) # m^3 m.rhow = Param(initialize=700.0) # kg/m^3 m.cpw = Param(initialize=3.1) # kJ/kg/K - m.Ca0 = Param(initialize=data['Ca0']) # kmol/m^3) - m.Cb0 = Param(initialize=data['Cb0']) # kmol/m^3) - m.Cc0 = Param(initialize=data['Cc0']) # kmol/m^3) + m.Ca0 = Param(initialize=data["Ca0"]) # kmol/m^3) + m.Cb0 = Param(initialize=data["Cb0"]) # kmol/m^3) + m.Cc0 = Param(initialize=data["Cc0"]) # kmol/m^3) m.Tr0 = Param(initialize=300.0) # K m.Vr0 = Param(initialize=1.0) # m^3 @@ -92,9 +101,9 @@ def generate_model(data): # def _initTc(m, t): if t < 10800: - return data['Tc1'] + return data["Tc1"] else: - return data['Tc2'] + return data["Tc2"] m.Tc = Param( m.time, initialize=_initTc, default=_initTc @@ -102,9 +111,9 @@ def _initTc(m, t): def _initFa(m, t): if t < 10800: - return data['Fa1'] + return data["Fa1"] else: - return data['Fa2'] + return data["Fa2"] m.Fa = Param( m.time, initialize=_initFa, default=_initFa @@ -230,7 +239,7 @@ def AllMeasurements(m): ) def MissingMeasurements(m): - if data['experiment'] == 1: + if data["experiment"] == 1: return sum( (m.Ca[t] - m.Ca_meas[t]) ** 2 + (m.Cb[t] - m.Cb_meas[t]) ** 2 @@ -238,7 +247,7 @@ def MissingMeasurements(m): + (m.Tr[t] - m.Tr_meas[t]) ** 2 for t in m.measT ) - elif data['experiment'] == 2: + elif data["experiment"] == 2: return sum((m.Tr[t] - m.Tr_meas[t]) ** 2 for t in m.measT) else: return sum( @@ -254,7 +263,7 @@ def total_cost_rule(model): m.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize) # Discretize model - disc = TransformationFactory('dae.collocation') + disc = TransformationFactory("dae.collocation") disc.apply_to(m, nfe=20, ncp=4) return m @@ -262,17 +271,17 @@ def total_cost_rule(model): def main(): # Data loaded from files file_dirname = dirname(abspath(str(__file__))) - file_name = abspath(join(file_dirname, 'exp2.out')) - with open(file_name, 'r') as infile: + file_name = abspath(join(file_dirname, "exp2.out")) + with open(file_name, "r") as infile: data = json.load(infile) - data['experiment'] = 2 + data["experiment"] = 2 model = generate_model(data) - solver = SolverFactory('ipopt') + solver = SolverFactory("ipopt") solver.solve(model) - print('k1 = ', model.k1()) - print('E1 = ', model.E1()) + print("k1 = ", model.k1()) + print("E1 = ", model.E1()) -if __name__ == '__main__': +if __name__ == "__main__": main() diff --git a/pyomo/contrib/parmest/graphics.py b/pyomo/contrib/parmest/graphics.py index f01622d2d17..b8dfa243b9a 100644 --- a/pyomo/contrib/parmest/graphics.py +++ b/pyomo/contrib/parmest/graphics.py @@ -29,7 +29,7 @@ # (e.g. python 3.5) get released that are either broken not # compatible, resulting in a SyntaxError sns, seaborn_available = attempt_import( - 'seaborn', catch_exceptions=(ImportError, SyntaxError) + "seaborn", catch_exceptions=(ImportError, SyntaxError) ) imports_available = ( @@ -93,17 +93,17 @@ def _get_data_slice(xvar, yvar, columns, data, theta_star): temp[col] = temp[col] + data[col].std() data = pd.concat([data, temp], ignore_index=True) - data_slice['obj'] = scipy.interpolate.griddata( + data_slice["obj"] = scipy.interpolate.griddata( np.array(data[columns]), - np.array(data[['obj']]), + np.array(data[["obj"]]), np.array(data_slice[columns]), - method='linear', + method="linear", rescale=True, ) X = data_slice[xvar] Y = data_slice[yvar] - Z = data_slice['obj'] + Z = data_slice["obj"] return X, Y, Z @@ -178,11 +178,11 @@ def _add_obj_contour(x, y, color, columns, data, theta_star, label=None): X, Y, Z = _get_data_slice(xvar, yvar, columns, data, theta_star) triang = matplotlib.tri.Triangulation(X, Y) - cmap = plt.cm.get_cmap('Greys') + cmap = matplotlib.colormaps["Greys"] plt.tricontourf(triang, Z, cmap=cmap) except: - print('Objective contour plot for', xvar, yvar, 'slice failed') + print("Objective contour plot for", xvar, yvar, "slice failed") def _set_axis_limits(g, axis_limits, theta_vals, theta_star): @@ -277,7 +277,7 @@ def pairwise_plot( assert isinstance(theta_star, (type(None), dict, pd.Series, pd.DataFrame)) assert isinstance(alpha, (type(None), int, float)) assert isinstance(distributions, list) - assert set(distributions).issubset(set(['MVN', 'KDE', 'Rect'])) + assert set(distributions).issubset(set(["MVN", "KDE", "Rect"])) assert isinstance(axis_limits, (type(None), dict)) assert isinstance(title, (type(None), str)) assert isinstance(add_obj_contour, bool) @@ -307,7 +307,7 @@ def pairwise_plot( theta_names = [ col for col in theta_values.columns - if (col not in ['obj']) + if (col not in ["obj"]) and (not isinstance(col, float)) and (not isinstance(col, int)) ] @@ -335,7 +335,7 @@ def pairwise_plot( g.map_diag(sns.distplot, kde=False, hist=True, norm_hist=False) # Plot filled contours using all theta values based on obj - if 'obj' in theta_values.columns and add_obj_contour: + if "obj" in theta_values.columns and add_obj_contour: g.map_offdiag( _add_obj_contour, columns=theta_names, @@ -349,10 +349,10 @@ def pairwise_plot( matplotlib.lines.Line2D( [0], [0], - marker='o', - color='w', - label='thetas', - markerfacecolor='cadetblue', + marker="o", + color="w", + label="thetas", + markerfacecolor="cadetblue", markersize=5, ) ) @@ -360,23 +360,23 @@ def pairwise_plot( # Plot theta* if theta_star is not None: g.map_offdiag( - _add_scatter, color='k', columns=theta_names, theta_star=theta_star + _add_scatter, color="k", columns=theta_names, theta_star=theta_star ) legend_elements.append( matplotlib.lines.Line2D( [0], [0], - marker='o', - color='w', - label='theta*', - markerfacecolor='k', + marker="o", + color="w", + label="theta*", + markerfacecolor="k", markersize=6, ) ) # Plot confidence regions - colors = ['r', 'mediumblue', 'darkgray'] + colors = ["r", "mediumblue", "darkgray"] if (alpha is not None) and (len(distributions) > 0): if theta_star is None: print( @@ -388,7 +388,7 @@ def pairwise_plot( mvn_dist = None kde_dist = None for i, dist in enumerate(distributions): - if dist == 'Rect': + if dist == "Rect": lb, ub = fit_rect_dist(thetas, alpha) g.map_offdiag( _add_rectangle_CI, @@ -401,7 +401,7 @@ def pairwise_plot( matplotlib.lines.Line2D([0], [0], color=colors[i], lw=1, label=dist) ) - elif dist == 'MVN': + elif dist == "MVN": mvn_dist = fit_mvn_dist(thetas) Z = mvn_dist.pdf(thetas) score = stats.scoreatpercentile(Z, (1 - alpha) * 100) @@ -418,7 +418,7 @@ def pairwise_plot( matplotlib.lines.Line2D([0], [0], color=colors[i], lw=1, label=dist) ) - elif dist == 'KDE': + elif dist == "KDE": kde_dist = fit_kde_dist(thetas) Z = kde_dist.pdf(thetas.transpose()) score = stats.scoreatpercentile(Z, (1 - alpha) * 100) @@ -438,12 +438,12 @@ def pairwise_plot( _set_axis_limits(g, axis_limits, thetas, theta_star) for ax in g.axes.flatten(): - ax.ticklabel_format(style='sci', scilimits=(-2, 2), axis='both') + ax.ticklabel_format(style="sci", scilimits=(-2, 2), axis="both") if add_legend: xvar, yvar, loc = _get_variables(ax, theta_names) if loc == (len(theta_names) - 1, 0): - ax.legend(handles=legend_elements, loc='best', prop={'size': 8}) + ax.legend(handles=legend_elements, loc="best", prop={"size": 8}) if title: g.fig.subplots_adjust(top=0.9) g.fig.suptitle(title) @@ -474,7 +474,7 @@ def pairwise_plot( ax.tick_params(reset=True) if add_legend: - ax.legend(handles=legend_elements, loc='best', prop={'size': 8}) + ax.legend(handles=legend_elements, loc="best", prop={"size": 8}) plt.close(g.fig) @@ -563,15 +563,15 @@ def _get_grouped_data(data1, data2, normalize, group_names): # Combine data1 and data2 to create a grouped histogram data = pd.concat({group_names[0]: data1, group_names[1]: data2}) data.reset_index(level=0, inplace=True) - data.rename(columns={'level_0': 'set'}, inplace=True) + data.rename(columns={"level_0": "set"}, inplace=True) - data = data.melt(id_vars='set', value_vars=data1.columns, var_name='columns') + data = data.melt(id_vars="set", value_vars=data1.columns, var_name="columns") return data def grouped_boxplot( - data1, data2, normalize=False, group_names=['data1', 'data2'], filename=None + data1, data2, normalize=False, group_names=["data1", "data2"], filename=None ): """ Plot a grouped boxplot to compare two datasets @@ -600,11 +600,11 @@ def grouped_boxplot( data = _get_grouped_data(data1, data2, normalize, group_names) plt.figure() - sns.boxplot(data=data, hue='set', y='value', x='columns', order=data1.columns) + sns.boxplot(data=data, hue="set", y="value", x="columns", order=data1.columns) - plt.gca().legend().set_title('') - plt.gca().set_xlabel('') - plt.gca().set_ylabel('') + plt.gca().legend().set_title("") + plt.gca().set_xlabel("") + plt.gca().set_ylabel("") if filename is None: plt.show() @@ -614,7 +614,7 @@ def grouped_boxplot( def grouped_violinplot( - data1, data2, normalize=False, group_names=['data1', 'data2'], filename=None + data1, data2, normalize=False, group_names=["data1", "data2"], filename=None ): """ Plot a grouped violinplot to compare two datasets @@ -644,12 +644,12 @@ def grouped_violinplot( plt.figure() sns.violinplot( - data=data, hue='set', y='value', x='columns', order=data1.columns, split=True + data=data, hue="set", y="value", x="columns", order=data1.columns, split=True ) - plt.gca().legend().set_title('') - plt.gca().set_xlabel('') - plt.gca().set_ylabel('') + plt.gca().legend().set_title("") + plt.gca().set_xlabel("") + plt.gca().set_ylabel("") if filename is None: plt.show() diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py index f26ecec2fce..2cc8ad36b0a 100644 --- a/pyomo/contrib/parmest/tests/test_parmest.py +++ b/pyomo/contrib/parmest/tests/test_parmest.py @@ -22,7 +22,7 @@ import platform -is_osx = platform.mac_ver()[0] != '' +is_osx = platform.mac_ver()[0] != "" import pyomo.common.unittest as unittest import sys @@ -38,11 +38,11 @@ from pyomo.opt import SolverFactory -ipopt_available = SolverFactory('ipopt').available() +ipopt_available = SolverFactory("ipopt").available() from pyomo.common.fileutils import find_library -pynumero_ASL_available = False if find_library('pynumero_ASL') is None else True +pynumero_ASL_available = False if find_library("pynumero_ASL") is None else True testdir = os.path.dirname(os.path.abspath(__file__)) @@ -61,10 +61,10 @@ def setUp(self): # Note, the data used in this test has been corrected to use data.loc[5,'hour'] = 7 (instead of 6) data = pd.DataFrame( data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]], - columns=['hour', 'y'], + columns=["hour", "y"], ) - theta_names = ['asymptote', 'rate_constant'] + theta_names = ["asymptote", "rate_constant"] def SSE(model, data): expr = sum( @@ -73,7 +73,7 @@ def SSE(model, data): ) return expr - solver_options = {'tol': 1e-8} + solver_options = {"tol": 1e-8} self.data = data self.pest = parmest.Estimator( @@ -90,10 +90,10 @@ def test_theta_est(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper @unittest.skipIf( @@ -105,14 +105,14 @@ def test_bootstrap(self): num_bootstraps = 10 theta_est = self.pest.theta_est_bootstrap(num_bootstraps, return_samples=True) - num_samples = theta_est['samples'].apply(len) + num_samples = theta_est["samples"].apply(len) self.assertTrue(len(theta_est.index), 10) self.assertTrue(num_samples.equals(pd.Series([6] * 10))) - del theta_est['samples'] + del theta_est["samples"] # apply confidence region test - CR = self.pest.confidence_region_test(theta_est, 'MVN', [0.5, 0.75, 1.0]) + CR = self.pest.confidence_region_test(theta_est, "MVN", [0.5, 0.75, 1.0]) self.assertTrue(set(CR.columns) >= set([0.5, 0.75, 1.0])) self.assertTrue(CR[0.5].sum() == 5) @@ -121,7 +121,7 @@ def test_bootstrap(self): graphics.pairwise_plot(theta_est) graphics.pairwise_plot(theta_est, thetavals) - graphics.pairwise_plot(theta_est, thetavals, 0.8, ['MVN', 'KDE', 'Rect']) + graphics.pairwise_plot(theta_est, thetavals, 0.8, ["MVN", "KDE", "Rect"]) @unittest.skipIf( not graphics.imports_available, "parmest.graphics imports are unavailable" @@ -151,7 +151,7 @@ def test_leaveNout(self): self.assertTrue(lNo_theta.shape == (6, 2)) results = self.pest.leaveNout_bootstrap_test( - 1, None, 3, 'Rect', [0.5, 1.0], seed=5436 + 1, None, 3, "Rect", [0.5, 1.0], seed=5436 ) self.assertTrue(len(results) == 6) # 6 lNo samples i = 1 @@ -221,10 +221,10 @@ def test_theta_est_cov(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper # Covariance matrix @@ -239,22 +239,22 @@ def test_theta_est_cov(self): ) # -0.4322 from paper self.assertAlmostEqual(cov.iloc[1, 1], 0.04124, places=2) # 0.04124 from paper - ''' Why does the covariance matrix from parmest not match the paper? Parmest is + """ Why does the covariance matrix from parmest not match the paper? Parmest is calculating the exact reduced Hessian. The paper (Rooney and Bielger, 2001) likely employed the first order approximation common for nonlinear regression. The paper values were verified with Scipy, which uses the same first order approximation. The formula used in parmest was verified against equations (7-5-15) and (7-5-16) in "Nonlinear Parameter Estimation", Y. Bard, 1974. - ''' + """ def test_cov_scipy_least_squares_comparison(self): - ''' + """ Scipy results differ in the 3rd decimal place from the paper. It is possible the paper used an alternative finite difference approximation for the Jacobian. - ''' + """ def model(theta, t): - ''' + """ Model to be fitted y = model(theta, t) Arguments: theta: vector of fitted parameters @@ -262,32 +262,32 @@ def model(theta, t): Returns: y: model predictions [need to check paper for units] - ''' + """ asymptote = theta[0] rate_constant = theta[1] return asymptote * (1 - np.exp(-rate_constant * t)) def residual(theta, t, y): - ''' + """ Calculate residuals Arguments: theta: vector of fitted parameters t: independent variable [hours] y: dependent variable [?] - ''' + """ return y - model(theta, t) # define data - t = self.data['hour'].to_numpy() - y = self.data['y'].to_numpy() + t = self.data["hour"].to_numpy() + y = self.data["y"].to_numpy() # define initial guess theta_guess = np.array([15, 0.5]) ## solve with optimize.least_squares sol = scipy.optimize.least_squares( - residual, theta_guess, method='trf', args=(t, y), verbose=2 + residual, theta_guess, method="trf", args=(t, y), verbose=2 ) theta_hat = sol.x @@ -313,18 +313,18 @@ def residual(theta, t, y): self.assertAlmostEqual(cov[1, 1], 0.04124, places=2) # 0.04124 from paper def test_cov_scipy_curve_fit_comparison(self): - ''' + """ Scipy results differ in the 3rd decimal place from the paper. It is possible the paper used an alternative finite difference approximation for the Jacobian. - ''' + """ ## solve with optimize.curve_fit def model(t, asymptote, rate_constant): return asymptote * (1 - np.exp(-rate_constant * t)) # define data - t = self.data['hour'].to_numpy() - y = self.data['y'].to_numpy() + t = self.data["hour"].to_numpy() + y = self.data["y"].to_numpy() # define initial guess theta_guess = np.array([15, 0.5]) @@ -351,7 +351,7 @@ class TestModelVariants(unittest.TestCase): def setUp(self): self.data = pd.DataFrame( data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]], - columns=['hour', 'y'], + columns=["hour", "y"], ) def rooney_biegler_params(data): @@ -371,16 +371,16 @@ def response_rule(m, h): def rooney_biegler_indexed_params(data): model = pyo.ConcreteModel() - model.param_names = pyo.Set(initialize=['asymptote', 'rate_constant']) + model.param_names = pyo.Set(initialize=["asymptote", "rate_constant"]) model.theta = pyo.Param( model.param_names, - initialize={'asymptote': 15, 'rate_constant': 0.5}, + initialize={"asymptote": 15, "rate_constant": 0.5}, mutable=True, ) def response_rule(m, h): - expr = m.theta['asymptote'] * ( - 1 - pyo.exp(-m.theta['rate_constant'] * h) + expr = m.theta["asymptote"] * ( + 1 - pyo.exp(-m.theta["rate_constant"] * h) ) return expr @@ -407,20 +407,20 @@ def response_rule(m, h): def rooney_biegler_indexed_vars(data): model = pyo.ConcreteModel() - model.var_names = pyo.Set(initialize=['asymptote', 'rate_constant']) + model.var_names = pyo.Set(initialize=["asymptote", "rate_constant"]) model.theta = pyo.Var( - model.var_names, initialize={'asymptote': 15, 'rate_constant': 0.5} + model.var_names, initialize={"asymptote": 15, "rate_constant": 0.5} ) model.theta[ - 'asymptote' + "asymptote" ].fixed = ( True # parmest will unfix theta variables, even when they are indexed ) - model.theta['rate_constant'].fixed = True + model.theta["rate_constant"].fixed = True def response_rule(m, h): - expr = m.theta['asymptote'] * ( - 1 - pyo.exp(-m.theta['rate_constant'] * h) + expr = m.theta["asymptote"] * ( + 1 - pyo.exp(-m.theta["rate_constant"] * h) ) return expr @@ -437,41 +437,41 @@ def SSE(model, data): self.objective_function = SSE - theta_vals = pd.DataFrame([20, 1], index=['asymptote', 'rate_constant']).T + theta_vals = pd.DataFrame([20, 1], index=["asymptote", "rate_constant"]).T theta_vals_index = pd.DataFrame( [20, 1], index=["theta['asymptote']", "theta['rate_constant']"] ).T self.input = { - 'param': { - 'model': rooney_biegler_params, - 'theta_names': ['asymptote', 'rate_constant'], - 'theta_vals': theta_vals, + "param": { + "model": rooney_biegler_params, + "theta_names": ["asymptote", "rate_constant"], + "theta_vals": theta_vals, }, - 'param_index': { - 'model': rooney_biegler_indexed_params, - 'theta_names': ['theta'], - 'theta_vals': theta_vals_index, + "param_index": { + "model": rooney_biegler_indexed_params, + "theta_names": ["theta"], + "theta_vals": theta_vals_index, }, - 'vars': { - 'model': rooney_biegler_vars, - 'theta_names': ['asymptote', 'rate_constant'], - 'theta_vals': theta_vals, + "vars": { + "model": rooney_biegler_vars, + "theta_names": ["asymptote", "rate_constant"], + "theta_vals": theta_vals, }, - 'vars_index': { - 'model': rooney_biegler_indexed_vars, - 'theta_names': ['theta'], - 'theta_vals': theta_vals_index, + "vars_index": { + "model": rooney_biegler_indexed_vars, + "theta_names": ["theta"], + "theta_vals": theta_vals_index, }, - 'vars_quoted_index': { - 'model': rooney_biegler_indexed_vars, - 'theta_names': ["theta['asymptote']", "theta['rate_constant']"], - 'theta_vals': theta_vals_index, + "vars_quoted_index": { + "model": rooney_biegler_indexed_vars, + "theta_names": ["theta['asymptote']", "theta['rate_constant']"], + "theta_vals": theta_vals_index, }, - 'vars_str_index': { - 'model': rooney_biegler_indexed_vars, - 'theta_names': ["theta[asymptote]", "theta[rate_constant]"], - 'theta_vals': theta_vals_index, + "vars_str_index": { + "model": rooney_biegler_indexed_vars, + "theta_names": ["theta[asymptote]", "theta[rate_constant]"], + "theta_vals": theta_vals_index, }, } @@ -483,9 +483,9 @@ def SSE(model, data): def test_parmest_basics(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) @@ -505,15 +505,15 @@ def test_parmest_basics(self): cov.iloc[1, 1], 0.04193591, places=2 ) # 0.04124 from paper - obj_at_theta = pest.objective_at_theta(parmest_input['theta_vals']) - self.assertAlmostEqual(obj_at_theta['obj'][0], 16.531953, places=2) + obj_at_theta = pest.objective_at_theta(parmest_input["theta_vals"]) + self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2) def test_parmest_basics_with_initialize_parmest_model_option(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) @@ -534,22 +534,22 @@ def test_parmest_basics_with_initialize_parmest_model_option(self): ) # 0.04124 from paper obj_at_theta = pest.objective_at_theta( - parmest_input['theta_vals'], initialize_parmest_model=True + parmest_input["theta_vals"], initialize_parmest_model=True ) - self.assertAlmostEqual(obj_at_theta['obj'][0], 16.531953, places=2) + self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2) def test_parmest_basics_with_square_problem_solve(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) obj_at_theta = pest.objective_at_theta( - parmest_input['theta_vals'], initialize_parmest_model=True + parmest_input["theta_vals"], initialize_parmest_model=True ) objval, thetavals, cov = pest.theta_est(calc_cov=True, cov_n=6) @@ -568,14 +568,14 @@ def test_parmest_basics_with_square_problem_solve(self): cov.iloc[1, 1], 0.04193591, places=2 ) # 0.04124 from paper - self.assertAlmostEqual(obj_at_theta['obj'][0], 16.531953, places=2) + self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2) def test_parmest_basics_with_square_problem_solve_no_theta_vals(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) @@ -632,17 +632,17 @@ def setUp(self): [1.90, 10000, 4491.3, 1049.4, 920.5, 1769.4], [1.95, 10000, 4538.8, 1045.8, 893.9, 1760.8], ], - columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd'], + columns=["sv", "caf", "ca", "cb", "cc", "cd"], ) - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] def SSE(model, data): expr = ( - (float(data['ca']) - model.ca) ** 2 - + (float(data['cb']) - model.cb) ** 2 - + (float(data['cc']) - model.cc) ** 2 - + (float(data['cd']) - model.cd) ** 2 + (float(data.iloc[0]["ca"]) - model.ca) ** 2 + + (float(data.iloc[0]["cb"]) - model.cb) ** 2 + + (float(data.iloc[0]["cc"]) - model.cc) ** 2 + + (float(data.iloc[0]["cd"]) - model.cd) ** 2 ) return expr @@ -656,13 +656,13 @@ def test_theta_est(self): # used in data reconciliation objval, thetavals = self.pest.theta_est() - self.assertAlmostEqual(thetavals['k1'], 5.0 / 6.0, places=4) - self.assertAlmostEqual(thetavals['k2'], 5.0 / 3.0, places=4) - self.assertAlmostEqual(thetavals['k3'], 1.0 / 6000.0, places=7) + self.assertAlmostEqual(thetavals["k1"], 5.0 / 6.0, places=4) + self.assertAlmostEqual(thetavals["k2"], 5.0 / 3.0, places=4) + self.assertAlmostEqual(thetavals["k3"], 1.0 / 6000.0, places=7) def test_return_values(self): objval, thetavals, data_rec = self.pest.theta_est( - return_values=['ca', 'cb', 'cc', 'cd', 'caf'] + return_values=["ca", "cb", "cc", "cd", "caf"] ) self.assertAlmostEqual(data_rec["cc"].loc[18], 893.84924, places=3) @@ -679,9 +679,9 @@ class TestReactorDesign_DAE(unittest.TestCase): def setUp(self): def ABC_model(data): - ca_meas = data['ca'] - cb_meas = data['cb'] - cc_meas = data['cc'] + ca_meas = data["ca"] + cb_meas = data["cb"] + cc_meas = data["cc"] if isinstance(data, pd.DataFrame): meas_t = data.index # time index @@ -754,7 +754,7 @@ def total_cost_rule(model): rule=total_cost_rule, sense=pyo.minimize ) - disc = pyo.TransformationFactory('dae.collocation') + disc = pyo.TransformationFactory("dae.collocation") disc.apply_to(m, nfe=20, ncp=2) return m @@ -785,15 +785,15 @@ def total_cost_rule(model): [4.737, 0.004, 0.036, 0.971], [5.000, -0.024, 0.028, 0.985], ] - data = pd.DataFrame(data, columns=['t', 'ca', 'cb', 'cc']) - data_df = data.set_index('t') + data = pd.DataFrame(data, columns=["t", "ca", "cb", "cc"]) + data_df = data.set_index("t") data_dict = { - 'ca': {k: v for (k, v) in zip(data.t, data.ca)}, - 'cb': {k: v for (k, v) in zip(data.t, data.cb)}, - 'cc': {k: v for (k, v) in zip(data.t, data.cc)}, + "ca": {k: v for (k, v) in zip(data.t, data.ca)}, + "cb": {k: v for (k, v) in zip(data.t, data.cb)}, + "cc": {k: v for (k, v) in zip(data.t, data.cc)}, } - theta_names = ['k1', 'k2'] + theta_names = ["k1", "k2"] self.pest_df = parmest.Estimator(ABC_model, [data_df], theta_names) self.pest_dict = parmest.Estimator(ABC_model, [data_dict], theta_names) @@ -815,30 +815,30 @@ def test_dataformats(self): obj2, theta2 = self.pest_dict.theta_est() self.assertAlmostEqual(obj1, obj2, places=6) - self.assertAlmostEqual(theta1['k1'], theta2['k1'], places=6) - self.assertAlmostEqual(theta1['k2'], theta2['k2'], places=6) + self.assertAlmostEqual(theta1["k1"], theta2["k1"], places=6) + self.assertAlmostEqual(theta1["k2"], theta2["k2"], places=6) def test_return_continuous_set(self): - ''' + """ test if ContinuousSet elements are returned correctly from theta_est() - ''' - obj1, theta1, return_vals1 = self.pest_df.theta_est(return_values=['time']) - obj2, theta2, return_vals2 = self.pest_dict.theta_est(return_values=['time']) - self.assertAlmostEqual(return_vals1['time'].loc[0][18], 2.368, places=3) - self.assertAlmostEqual(return_vals2['time'].loc[0][18], 2.368, places=3) + """ + obj1, theta1, return_vals1 = self.pest_df.theta_est(return_values=["time"]) + obj2, theta2, return_vals2 = self.pest_dict.theta_est(return_values=["time"]) + self.assertAlmostEqual(return_vals1["time"].loc[0][18], 2.368, places=3) + self.assertAlmostEqual(return_vals2["time"].loc[0][18], 2.368, places=3) def test_return_continuous_set_multiple_datasets(self): - ''' + """ test if ContinuousSet elements are returned correctly from theta_est() - ''' + """ obj1, theta1, return_vals1 = self.pest_df_multiple.theta_est( - return_values=['time'] + return_values=["time"] ) obj2, theta2, return_vals2 = self.pest_dict_multiple.theta_est( - return_values=['time'] + return_values=["time"] ) - self.assertAlmostEqual(return_vals1['time'].loc[1][18], 2.368, places=3) - self.assertAlmostEqual(return_vals2['time'].loc[1][18], 2.368, places=3) + self.assertAlmostEqual(return_vals1["time"].loc[1][18], 2.368, places=3) + self.assertAlmostEqual(return_vals2["time"].loc[1][18], 2.368, places=3) def test_covariance(self): from pyomo.contrib.interior_point.inverse_reduced_hessian import ( @@ -862,13 +862,13 @@ def test_covariance(self): l = len(vars_list) cov_interior_point = 2 * obj / (n - l) * inv_red_hes cov_interior_point = pd.DataFrame( - cov_interior_point, ['k1', 'k2'], ['k1', 'k2'] + cov_interior_point, ["k1", "k2"], ["k1", "k2"] ) cov_diff = (cov - cov_interior_point).abs().sum().sum() - self.assertTrue(cov.loc['k1', 'k1'] > 0) - self.assertTrue(cov.loc['k2', 'k2'] > 0) + self.assertTrue(cov.loc["k1", "k1"] > 0) + self.assertTrue(cov.loc["k2", "k2"] > 0) self.assertAlmostEqual(cov_diff, 0, places=6) @@ -886,10 +886,10 @@ def setUp(self): # Note, the data used in this test has been corrected to use data.loc[5,'hour'] = 7 (instead of 6) data = pd.DataFrame( data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]], - columns=['hour', 'y'], + columns=["hour", "y"], ) - theta_names = ['asymptote', 'rate_constant'] + theta_names = ["asymptote", "rate_constant"] def SSE(model, data): expr = sum( @@ -898,7 +898,7 @@ def SSE(model, data): ) return expr - solver_options = {'tol': 1e-8} + solver_options = {"tol": 1e-8} self.data = data self.pest = parmest.Estimator( @@ -916,15 +916,15 @@ def test_theta_est_with_square_initialization(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper def test_theta_est_with_square_initialization_and_custom_init_theta(self): theta_vals_init = pd.DataFrame( - data=[[19.0, 0.5]], columns=['asymptote', 'rate_constant'] + data=[[19.0, 0.5]], columns=["asymptote", "rate_constant"] ) obj_init = self.pest.objective_at_theta( theta_values=theta_vals_init, initialize_parmest_model=True @@ -932,10 +932,10 @@ def test_theta_est_with_square_initialization_and_custom_init_theta(self): objval, thetavals = self.pest.theta_est() self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper def test_theta_est_with_square_initialization_diagnostic_mode_true(self): @@ -945,14 +945,14 @@ def test_theta_est_with_square_initialization_diagnostic_mode_true(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper self.pest.diagnostic_mode = False -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/parmest/tests/test_scenariocreator.py b/pyomo/contrib/parmest/tests/test_scenariocreator.py index a2dcf4c2739..22a851ae32e 100644 --- a/pyomo/contrib/parmest/tests/test_scenariocreator.py +++ b/pyomo/contrib/parmest/tests/test_scenariocreator.py @@ -24,7 +24,7 @@ import pyomo.environ as pyo from pyomo.environ import SolverFactory -ipopt_available = SolverFactory('ipopt').available() +ipopt_available = SolverFactory("ipopt").available() testdir = os.path.dirname(os.path.abspath(__file__)) @@ -63,17 +63,17 @@ def setUp(self): [1.90, 10000, 4491.3, 1049.4, 920.5, 1769.4], [1.95, 10000, 4538.8, 1045.8, 893.9, 1760.8], ], - columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd'], + columns=["sv", "caf", "ca", "cb", "cc", "cd"], ) - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] def SSE(model, data): expr = ( - (float(data['ca']) - model.ca) ** 2 - + (float(data['cb']) - model.cb) ** 2 - + (float(data['cc']) - model.cc) ** 2 - + (float(data['cd']) - model.cd) ** 2 + (float(data.iloc[0]["ca"]) - model.ca) ** 2 + + (float(data.iloc[0]["cb"]) - model.cb) ** 2 + + (float(data.iloc[0]["cc"]) - model.cc) ** 2 + + (float(data.iloc[0]["cd"]) - model.cd) ** 2 ) return expr @@ -116,7 +116,7 @@ def setUp(self): import json # Vars to estimate in parmest - theta_names = ['k1', 'k2', 'E1', 'E2'] + theta_names = ["k1", "k2", "E1", "E2"] self.fbase = os.path.join(testdir, "..", "examples", "semibatch") # Data, list of dictionaries @@ -124,7 +124,7 @@ def setUp(self): for exp_num in range(10): fname = "exp" + str(exp_num + 1) + ".out" fullname = os.path.join(self.fbase, fname) - with open(fullname, 'r') as infile: + with open(fullname, "r") as infile: d = json.load(infile) data.append(d) @@ -142,5 +142,5 @@ def test_semibatch_bootstrap(self): self.assertAlmostEqual(tval, 20.64, places=1) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/parmest/tests/test_utils.py b/pyomo/contrib/parmest/tests/test_utils.py index bd0706ac38d..514c14b1e82 100644 --- a/pyomo/contrib/parmest/tests/test_utils.py +++ b/pyomo/contrib/parmest/tests/test_utils.py @@ -16,7 +16,7 @@ import pyomo.contrib.parmest.parmest as parmest from pyomo.opt import SolverFactory -ipopt_available = SolverFactory('ipopt').available() +ipopt_available = SolverFactory("ipopt").available() @unittest.skipIf( @@ -45,13 +45,13 @@ def test_convert_param_to_var(self): [1.10, 10000, 3535.1, 1064.8, 1613.3, 1893.4], [1.15, 10000, 3609.1, 1067.8, 1547.5, 1887.8], ], - columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd'], + columns=["sv", "caf", "ca", "cb", "cc", "cd"], ) - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] instance = reactor_design_model(data.loc[0]) - solver = pyo.SolverFactory('ipopt') + solver = pyo.SolverFactory("ipopt") solver.solve(instance) instance_vars = parmest.utils.convert_params_to_vars( diff --git a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py index 766ef96322a..cedbf430a12 100644 --- a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py +++ b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py @@ -289,6 +289,13 @@ class PyomoCyIpoptSolver(object): description="Set the function that will be called each iteration.", ), ) + CONFIG.declare( + "halt_on_evaluation_error", + ConfigValue( + default=None, + description="Whether to halt if a function or derivative evaluation fails", + ), + ) def __init__(self, **kwds): """Create an instance of the CyIpoptSolver. You must @@ -332,7 +339,9 @@ def solve(self, model, **kwds): nlp = pyomo_nlp.PyomoNLP(model) problem = cyipopt_interface.CyIpoptNLP( - nlp, intermediate_callback=config.intermediate_callback + nlp, + intermediate_callback=config.intermediate_callback, + halt_on_evaluation_error=config.halt_on_evaluation_error, ) ng = len(problem.g_lb()) nx = len(problem.x_lb()) diff --git a/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py b/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py index 2a7edb430d4..7ead30117cb 100644 --- a/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py +++ b/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py @@ -24,6 +24,7 @@ if not (numpy_available and scipy_available): raise unittest.SkipTest("Pynumero needs scipy and numpy to run NLP tests") +from pyomo.contrib.pynumero.exceptions import PyNumeroEvaluationError from pyomo.contrib.pynumero.asl import AmplInterface if not AmplInterface.available(): @@ -34,12 +35,18 @@ from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.interfaces.cyipopt_interface import ( + cyipopt, cyipopt_available, CyIpoptNLP, ) from pyomo.contrib.pynumero.algorithms.solvers.cyipopt_solver import CyIpoptSolver +if cyipopt_available: + # We don't raise unittest.SkipTest if not cyipopt_available as there is a + # test below that tests an exception when cyipopt is unavailable. + cyipopt_ge_1_3 = hasattr(cyipopt, "CyIpoptEvaluationError") + def create_model1(): m = pyo.ConcreteModel() @@ -155,6 +162,29 @@ def f(model): return model +def make_hs071_model(): + # This is a model that is mathematically equivalent to the Hock-Schittkowski + # test problem 071, but that will trigger an evaluation error if x[0] goes + # above 1.1. + m = pyo.ConcreteModel() + m.x = pyo.Var([0, 1, 2, 3], bounds=(1.0, 5.0)) + m.x[0] = 1.0 + m.x[1] = 5.0 + m.x[2] = 5.0 + m.x[3] = 1.0 + m.obj = pyo.Objective(expr=m.x[0] * m.x[3] * (m.x[0] + m.x[1] + m.x[2]) + m.x[2]) + # This expression evaluates to zero, but is not well defined when x[0] > 1.1 + trivial_expr_with_eval_error = (pyo.sqrt(1.1 - m.x[0])) ** 2 + m.x[0] - 1.1 + m.ineq1 = pyo.Constraint(expr=m.x[0] * m.x[1] * m.x[2] * m.x[3] >= 25.0) + m.eq1 = pyo.Constraint( + expr=( + m.x[0] ** 2 + m.x[1] ** 2 + m.x[2] ** 2 + m.x[3] ** 2 + == 40.0 + trivial_expr_with_eval_error + ) + ) + return m + + @unittest.skipIf(cyipopt_available, "cyipopt is available") class TestCyIpoptNotAvailable(unittest.TestCase): def test_not_available_exception(self): @@ -257,3 +287,32 @@ def test_options(self): x, info = solver.solve(tee=False) nlp.set_primals(x) self.assertAlmostEqual(nlp.evaluate_objective(), -5.0879028e02, places=5) + + @unittest.skipUnless( + cyipopt_available and cyipopt_ge_1_3, "cyipopt version < 1.3.0" + ) + def test_hs071_evalerror(self): + m = make_hs071_model() + solver = pyo.SolverFactory("cyipopt") + res = solver.solve(m, tee=True) + + x = list(m.x[:].value) + expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829]) + np.testing.assert_allclose(x, expected_x) + + def test_hs071_evalerror_halt(self): + m = make_hs071_model() + solver = pyo.SolverFactory("cyipopt", halt_on_evaluation_error=True) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + res = solver.solve(m, tee=True) + + @unittest.skipIf( + not cyipopt_available or cyipopt_ge_1_3, "cyipopt version >= 1.3.0" + ) + def test_hs071_evalerror_old_cyipopt(self): + m = make_hs071_model() + solver = pyo.SolverFactory("cyipopt") + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + res = solver.solve(m, tee=True) diff --git a/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py b/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py index 19e74625d03..fc9c45c6d1a 100644 --- a/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py +++ b/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py @@ -23,6 +23,7 @@ import abc from pyomo.common.dependencies import attempt_import, numpy as np, numpy_available +from pyomo.contrib.pynumero.exceptions import PyNumeroEvaluationError def _cyipopt_importer(): @@ -252,7 +253,7 @@ def intermediate( class CyIpoptNLP(CyIpoptProblemInterface): - def __init__(self, nlp, intermediate_callback=None): + def __init__(self, nlp, intermediate_callback=None, halt_on_evaluation_error=None): """This class provides a CyIpoptProblemInterface for use with the CyIpoptSolver class that can take in an NLP as long as it provides vectors as numpy ndarrays and @@ -263,6 +264,23 @@ def __init__(self, nlp, intermediate_callback=None): self._nlp = nlp self._intermediate_callback = intermediate_callback + cyipopt_has_eval_error = cyipopt_available and hasattr( + cyipopt, "CyIpoptEvaluationError" + ) + if halt_on_evaluation_error is None: + # If using cyipopt >= 1.3, the default is to continue. + # Otherwise, the default is to halt (because we are forced to). + # + # If CyIpopt is not available, we "halt" (re-raise the original + # exception). + self._halt_on_evaluation_error = not cyipopt_has_eval_error + elif not halt_on_evaluation_error and not cyipopt_has_eval_error: + raise ValueError( + "halt_on_evaluation_error=False is only supported for cyipopt >= 1.3.0" + ) + else: + self._halt_on_evaluation_error = halt_on_evaluation_error + x = nlp.init_primals() y = nlp.init_duals() if np.any(np.isnan(y)): @@ -328,24 +346,54 @@ def scaling_factors(self): return obj_scaling, x_scaling, g_scaling def objective(self, x): - self._set_primals_if_necessary(x) - return self._nlp.evaluate_objective() + try: + self._set_primals_if_necessary(x) + return self._nlp.evaluate_objective() + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in objective function evaluation" + ) def gradient(self, x): - self._set_primals_if_necessary(x) - return self._nlp.evaluate_grad_objective() + try: + self._set_primals_if_necessary(x) + return self._nlp.evaluate_grad_objective() + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in objective gradient evaluation" + ) def constraints(self, x): - self._set_primals_if_necessary(x) - return self._nlp.evaluate_constraints() + try: + self._set_primals_if_necessary(x) + return self._nlp.evaluate_constraints() + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError("Error in constraint evaluation") def jacobianstructure(self): return self._jac_g.row, self._jac_g.col def jacobian(self, x): - self._set_primals_if_necessary(x) - self._nlp.evaluate_jacobian(out=self._jac_g) - return self._jac_g.data + try: + self._set_primals_if_necessary(x) + self._nlp.evaluate_jacobian(out=self._jac_g) + return self._jac_g.data + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in constraint Jacobian evaluation" + ) def hessianstructure(self): if not self._hessian_available: @@ -359,12 +407,20 @@ def hessian(self, x, y, obj_factor): if not self._hessian_available: raise ValueError("Hessian requested, but not supported by the NLP") - self._set_primals_if_necessary(x) - self._set_duals_if_necessary(y) - self._set_obj_factor_if_necessary(obj_factor) - self._nlp.evaluate_hessian_lag(out=self._hess_lag) - data = np.compress(self._hess_lower_mask, self._hess_lag.data) - return data + try: + self._set_primals_if_necessary(x) + self._set_duals_if_necessary(y) + self._set_obj_factor_if_necessary(obj_factor) + self._nlp.evaluate_hessian_lag(out=self._hess_lag) + data = np.compress(self._hess_lower_mask, self._hess_lag.data) + return data + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in Lagrangian Hessian evaluation" + ) def intermediate( self, diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py b/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py index 2c5d8ff7e4e..f28b7b9b549 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py +++ b/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py @@ -10,6 +10,7 @@ # ___________________________________________________________________________ import pyomo.common.unittest as unittest +import pyomo.environ as pyo from pyomo.contrib.pynumero.dependencies import ( numpy as np, @@ -20,20 +21,27 @@ if not (numpy_available and scipy_available): raise unittest.SkipTest("Pynumero needs scipy and numpy to run CyIpopt tests") +from pyomo.contrib.pynumero.exceptions import PyNumeroEvaluationError from pyomo.contrib.pynumero.asl import AmplInterface if not AmplInterface.available(): raise unittest.SkipTest("Pynumero needs the ASL extension to run CyIpopt tests") +from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.interfaces.cyipopt_interface import ( + cyipopt, cyipopt_available, CyIpoptProblemInterface, + CyIpoptNLP, ) if not cyipopt_available: raise unittest.SkipTest("CyIpopt is not available") +cyipopt_ge_1_3 = hasattr(cyipopt, "CyIpoptEvaluationError") + + class TestSubclassCyIpoptInterface(unittest.TestCase): def test_subclass_no_init(self): class MyCyIpoptProblem(CyIpoptProblemInterface): @@ -88,5 +96,129 @@ def hessian(self, x, y, obj_factor): problem.solve(x0) +def _get_model_nlp_interface(halt_on_evaluation_error=None): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], initialize=1.0) + m.obj = pyo.Objective(expr=m.x[1] * pyo.sqrt(m.x[2]) + m.x[1] * m.x[3]) + m.eq1 = pyo.Constraint(expr=m.x[1] * pyo.sqrt(m.x[2]) == 1.0) + nlp = PyomoNLP(m) + interface = CyIpoptNLP(nlp, halt_on_evaluation_error=halt_on_evaluation_error) + bad_primals = np.array([1.0, -2.0, 3.0]) + indices = nlp.get_primal_indices([m.x[1], m.x[2], m.x[3]]) + bad_primals = bad_primals[indices] + return m, nlp, interface, bad_primals + + +class TestCyIpoptVersionDependentConfig(unittest.TestCase): + @unittest.skipIf(cyipopt_ge_1_3, "cyipopt version >= 1.3.0") + def test_config_error(self): + _, nlp, _, _ = _get_model_nlp_interface() + with self.assertRaisesRegex(ValueError, "halt_on_evaluation_error"): + interface = CyIpoptNLP(nlp, halt_on_evaluation_error=False) + + @unittest.skipIf(cyipopt_ge_1_3, "cyipopt version >= 1.3.0") + def test_default_config_with_old_cyipopt(self): + _, nlp, _, bad_x = _get_model_nlp_interface() + interface = CyIpoptNLP(nlp) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.objective(bad_x) + + @unittest.skipIf(not cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_default_config_with_new_cyipopt(self): + _, nlp, _, bad_x = _get_model_nlp_interface() + interface = CyIpoptNLP(nlp) + msg = "Error in objective function" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.objective(bad_x) + + +class TestCyIpoptEvaluationErrors(unittest.TestCase): + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_objective(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in objective function" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.objective(bad_x) + + def test_error_in_objective_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.objective(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_gradient(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in objective gradient" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.gradient(bad_x) + + def test_error_in_gradient_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.gradient(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_constraints(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in constraint evaluation" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.constraints(bad_x) + + def test_error_in_constraints_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.constraints(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_jacobian(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in constraint Jacobian" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.jacobian(bad_x) + + def test_error_in_jacobian_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.jacobian(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_hessian(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in Lagrangian Hessian" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.hessian(bad_x, [1.0], 0.0) + + def test_error_in_hessian_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.hessian(bad_x, [1.0], 0.0) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/pyros/CHANGELOG.txt b/pyomo/contrib/pyros/CHANGELOG.txt index b1866ed955c..7977c37fb95 100644 --- a/pyomo/contrib/pyros/CHANGELOG.txt +++ b/pyomo/contrib/pyros/CHANGELOG.txt @@ -2,6 +2,18 @@ PyROS CHANGELOG =============== +------------------------------------------------------------------------------- +PyROS 1.2.8 12 Oct 2023 +------------------------------------------------------------------------------- +- Refactor PyROS separation routine, fix scenario selection heuristic +- Add efficiency for discrete uncertainty set separation +- Fix coefficient matching routine +- Fix subproblem timers and time accumulators +- Update and document PyROS solver logging system +- Fix iteration overcounting in event of `max_iter` termination status +- Fixes to (assembly of) PyROS `ROSolveResults` object + + ------------------------------------------------------------------------------- PyROS 1.2.7 26 Apr 2023 ------------------------------------------------------------------------------- diff --git a/pyomo/contrib/pyros/master_problem_methods.py b/pyomo/contrib/pyros/master_problem_methods.py index a0e2245cab1..dc4b6b957bb 100644 --- a/pyomo/contrib/pyros/master_problem_methods.py +++ b/pyomo/contrib/pyros/master_problem_methods.py @@ -22,7 +22,6 @@ adjust_solver_time_settings, revert_solver_max_time_adjustment, get_main_elapsed_time, - output_logger, ) from pyomo.contrib.pyros.solve_data import MasterProblemData, MasterResult from pyomo.opt.results import check_optimal_termination @@ -236,6 +235,13 @@ def solve_master_feasibility_problem(model_data, config): """ model = construct_master_feasibility_problem(model_data, config) + active_obj = next(model.component_data_objects(Objective, active=True)) + + config.progress_logger.debug("Solving master feasibility problem") + config.progress_logger.debug( + f" Initial objective (total slack): {value(active_obj)}" + ) + if config.solve_master_globally: solver = config.global_solver else: @@ -245,6 +251,7 @@ def solve_master_feasibility_problem(model_data, config): orig_setting, custom_setting_present = adjust_solver_time_settings( model_data.timing, solver, config ) + model_data.timing.start_timer("main.master_feasibility") timer.tic(msg=None) try: results = solver.solve(model, tee=config.tee, load_solutions=False) @@ -253,13 +260,14 @@ def solve_master_feasibility_problem(model_data, config): # (such as segmentation faults, function evaluation # errors, etc.) config.progress_logger.error( - f"Solver {repr(solver)} encountered exception attempting to " - "optimize master feasibility problem in iteration " - f"{model_data.iteration}" + f"Optimizer {repr(solver)} encountered exception " + "attempting to solve master feasibility problem in iteration " + f"{model_data.iteration}." ) raise else: setattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR, timer.toc(msg=None)) + model_data.timing.stop_timer("main.master_feasibility") finally: revert_solver_max_time_adjustment( solver, orig_setting, custom_setting_present, config @@ -273,6 +281,24 @@ def solve_master_feasibility_problem(model_data, config): } if results.solver.termination_condition in feasible_terminations: model.solutions.load_from(results) + config.progress_logger.debug( + f" Final objective (total slack): {value(active_obj)}" + ) + config.progress_logger.debug( + f" Termination condition: {results.solver.termination_condition}" + ) + config.progress_logger.debug( + f" Solve time: {getattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR)}s" + ) + else: + config.progress_logger.warning( + "Could not successfully solve master feasibility problem " + f"of iteration {model_data.iteration} with primary subordinate " + f"{'global' if config.solve_master_globally else 'local'} solver " + "to acceptable level. " + f"Termination stats:\n{results.solver}\n" + "Maintaining unoptimized point for master problem initialization." + ) # load master feasibility point to master model for master_var, feas_var in model_data.feasibility_problem_varmap: @@ -298,6 +324,9 @@ def minimize_dr_vars(model_data, config): ------- results : SolverResults Subordinate solver results for the polishing problem. + polishing_successful : bool + True if polishing model was solved to acceptable level, + False otherwise. """ # config.progress_logger.info("Executing decision rule variable polishing solve.") model = model_data.master_model @@ -467,11 +496,21 @@ def minimize_dr_vars(model_data, config): else: solver = config.local_solver + config.progress_logger.debug("Solving DR polishing problem") + + # NOTE: this objective evalaution may not be accurate, due + # to the current initialization scheme for the auxiliary + # variables. new initialization will be implemented in the + # near future. + polishing_obj = polishing_model.scenarios[0, 0].polishing_obj + config.progress_logger.debug(f" Initial DR norm: {value(polishing_obj)}") + # === Solve the polishing model timer = TicTocTimer() orig_setting, custom_setting_present = adjust_solver_time_settings( model_data.timing, solver, config ) + model_data.timing.start_timer("main.dr_polishing") timer.tic(msg=None) try: results = solver.solve(polishing_model, tee=config.tee, load_solutions=False) @@ -484,16 +523,35 @@ def minimize_dr_vars(model_data, config): raise else: setattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR, timer.toc(msg=None)) + model_data.timing.stop_timer("main.dr_polishing") finally: revert_solver_max_time_adjustment( solver, orig_setting, custom_setting_present, config ) + # interested in the time and termination status for debugging + # purposes + config.progress_logger.debug(" Done solving DR polishing problem") + config.progress_logger.debug( + f" Solve time: {getattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR)} s" + ) + config.progress_logger.debug( + f" Termination status: {results.solver.termination_condition} " + ) + # === Process solution by termination condition acceptable = {tc.globallyOptimal, tc.optimal, tc.locallyOptimal, tc.feasible} if results.solver.termination_condition not in acceptable: # continue with "unpolished" master model solution - return results + config.progress_logger.warning( + "Could not successfully solve DR polishing problem " + f"of iteration {model_data.iteration} with primary subordinate " + f"{'global' if config.solve_master_globally else 'local'} solver " + "to acceptable level. " + f"Termination stats:\n{results.solver}\n" + "Maintaining unpolished master problem solution." + ) + return results, False # update master model second-stage, state, and decision rule # variables to polishing model solution @@ -520,7 +578,34 @@ def minimize_dr_vars(model_data, config): for mvar, pvar in zip(master_dr.values(), polish_dr.values()): mvar.set_value(value(pvar), skip_validation=True) - return results + config.progress_logger.debug(f" Optimized DR norm: {value(polishing_obj)}") + config.progress_logger.debug(" Polished Master objective:") + + # print master solution + if config.objective_focus == ObjectiveType.worst_case: + worst_blk_idx = max( + model_data.master_model.scenarios.keys(), + key=lambda idx: value( + model_data.master_model.scenarios[idx].second_stage_objective + ), + ) + else: + worst_blk_idx = (0, 0) + + # debugging: summarize objective breakdown + worst_master_blk = model_data.master_model.scenarios[worst_blk_idx] + config.progress_logger.debug( + " First-stage objective: " f"{value(worst_master_blk.first_stage_objective)}" + ) + config.progress_logger.debug( + " Second-stage objective: " f"{value(worst_master_blk.second_stage_objective)}" + ) + polished_master_obj = value( + worst_master_blk.first_stage_objective + worst_master_blk.second_stage_objective + ) + config.progress_logger.debug(f" Objective: {polished_master_obj}") + + return results, True def add_p_robust_constraint(model_data, config): @@ -628,18 +713,27 @@ def solver_call_master(model_data, config, solver, solve_data): solver_term_cond_dict = {} if config.solve_master_globally: - backup_solvers = deepcopy(config.backup_global_solvers) + solvers = [solver] + config.backup_global_solvers else: - backup_solvers = deepcopy(config.backup_local_solvers) - backup_solvers.insert(0, solver) + solvers = [solver] + config.backup_local_solvers higher_order_decision_rule_efficiency(config, model_data) + solve_mode = "global" if config.solve_master_globally else "local" + config.progress_logger.debug("Solving master problem") + timer = TicTocTimer() - for opt in backup_solvers: + for idx, opt in enumerate(solvers): + if idx > 0: + config.progress_logger.warning( + f"Invoking backup solver {opt!r} " + f"(solver {idx + 1} of {len(solvers)}) for " + f"master problem of iteration {model_data.iteration}." + ) orig_setting, custom_setting_present = adjust_solver_time_settings( model_data.timing, opt, config ) + model_data.timing.start_timer("main.master") timer.tic(msg=None) try: results = opt.solve( @@ -653,12 +747,14 @@ def solver_call_master(model_data, config, solver, solve_data): # (such as segmentation faults, function evaluation # errors, etc.) config.progress_logger.error( - f"Solver {repr(opt)} encountered exception attempting to " - f"optimize master problem in iteration {model_data.iteration}" + f"Optimizer {repr(opt)} ({idx + 1} of {len(solvers)}) " + "encountered exception attempting to " + f"solve master problem in iteration {model_data.iteration}" ) raise else: setattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR, timer.toc(msg=None)) + model_data.timing.stop_timer("main.master") finally: revert_solver_max_time_adjustment( solver, orig_setting, custom_setting_present, config @@ -715,6 +811,25 @@ def solver_call_master(model_data, config, solver, solve_data): nlp_model.scenarios[0, 0].first_stage_objective ) + # debugging: log breakdown of master objective + config.progress_logger.debug(" Optimized master objective breakdown:") + config.progress_logger.debug( + f" First-stage objective: {master_soln.first_stage_objective}" + ) + config.progress_logger.debug( + f" Second-stage objective: {master_soln.second_stage_objective}" + ) + master_obj = ( + master_soln.first_stage_objective + master_soln.second_stage_objective + ) + config.progress_logger.debug(f" Objective: {master_obj}") + config.progress_logger.debug( + f" Termination condition: {results.solver.termination_condition}" + ) + config.progress_logger.debug( + f" Solve time: {getattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR)}s" + ) + master_soln.nominal_block = nlp_model.scenarios[0, 0] master_soln.results = results master_soln.master_model = nlp_model @@ -731,7 +846,6 @@ def solver_call_master(model_data, config, solver, solve_data): master_soln.pyros_termination_condition = ( pyrosTerminationCondition.time_out ) - output_logger(config=config, time_out=True, elapsed=elapsed) if not try_backup: return master_soln @@ -742,8 +856,9 @@ def solver_call_master(model_data, config, solver, solve_data): # NOTE: subproblem is written with variables set to their # initial values (not the final subsolver iterate) save_dir = config.subproblem_file_directory + serialization_msg = "" if save_dir and config.keepfiles: - name = os.path.join( + output_problem_path = os.path.join( save_dir, ( config.uncertainty_set.type @@ -754,15 +869,37 @@ def solver_call_master(model_data, config, solver, solve_data): + ".bar" ), ) - nlp_model.write(name, io_options={'symbolic_solver_labels': True}) - output_logger( - config=config, - master_error=True, - status_dict=solver_term_cond_dict, - filename=name, - iteration=model_data.iteration, + nlp_model.write( + output_problem_path, io_options={'symbolic_solver_labels': True} ) + serialization_msg = ( + " For debugging, problem has been serialized to the file " + f"{output_problem_path!r}." + ) + + deterministic_model_qual = ( + " (i.e., the deterministic model)" if model_data.iteration == 0 else "" + ) + deterministic_msg = ( + ( + " Please ensure your deterministic model " + f"is solvable by at least one of the subordinate {solve_mode} " + "optimizers provided." + ) + if model_data.iteration == 0 + else "" + ) master_soln.pyros_termination_condition = pyrosTerminationCondition.subsolver_error + config.progress_logger.warning( + f"Could not successfully solve master problem of iteration " + f"{model_data.iteration}{deterministic_model_qual} with any of the " + f"provided subordinate {solve_mode} optimizers. " + f"(Termination statuses: " + f"{[term_cond for term_cond in solver_term_cond_dict.values()]}.)" + f"{deterministic_msg}" + f"{serialization_msg}" + ) + return master_soln @@ -798,10 +935,6 @@ def solve_master(model_data, config): None, pyrosTerminationCondition.time_out, ) - - # log time out message - output_logger(config=config, time_out=True, elapsed=elapsed) - return master_soln solver = ( diff --git a/pyomo/contrib/pyros/pyros.py b/pyomo/contrib/pyros/pyros.py index 34db54b64e6..5b37b114722 100644 --- a/pyomo/contrib/pyros/pyros.py +++ b/pyomo/contrib/pyros/pyros.py @@ -37,15 +37,52 @@ transform_to_standard_form, turn_bounds_to_constraints, replace_uncertain_bounds_with_constraints, - output_logger, + IterationLogRecord, + setup_pyros_logger, + TimingData, ) from pyomo.contrib.pyros.solve_data import ROSolveResults from pyomo.contrib.pyros.pyros_algorithm_methods import ROSolver_iterative_solve from pyomo.contrib.pyros.uncertainty_sets import uncertainty_sets from pyomo.core.base import Constraint +from datetime import datetime -__version__ = "1.2.7" + +__version__ = "1.2.8" + + +default_pyros_solver_logger = setup_pyros_logger() + + +def _get_pyomo_version_info(): + """ + Get Pyomo version information. + """ + import os + import subprocess + from pyomo.version import version + + pyomo_version = version + commit_hash = "unknown" + + pyros_dir = os.path.join(*os.path.split(__file__)[:-1]) + commit_hash_command_args = [ + "git", + "-C", + f"{pyros_dir}", + "rev-parse", + "--short", + "HEAD", + ] + try: + commit_hash = ( + subprocess.check_output(commit_hash_command_args).decode("ascii").strip() + ) + except subprocess.CalledProcessError: + commit_hash = "unknown" + + return {"Pyomo version": pyomo_version, "Commit hash": commit_hash} def NonNegIntOrMinusOne(obj): @@ -485,13 +522,16 @@ def pyros_config(): CONFIG.declare( "progress_logger", PyROSConfigValue( - default="pyomo.contrib.pyros", + default=default_pyros_solver_logger, domain=a_logger, doc=( """ Logger (or name thereof) used for reporting PyROS solver - progress. If a `str` is specified, then - ``logging.getLogger(progress_logger)`` is used. + progress. If a `str` is specified, then ``progress_logger`` + is cast to ``logging.getLogger(progress_logger)``. + In the default case, `progress_logger` is set to + a :class:`pyomo.contrib.pyros.util.PreformattedLogger` + object of level ``logging.INFO``. """ ), is_optional=True, @@ -642,6 +682,7 @@ class PyROS(object): ''' CONFIG = pyros_config() + _LOG_LINE_LENGTH = 78 def available(self, exception_flag=True): """Check if solver is available.""" @@ -663,6 +704,138 @@ def __enter__(self): def __exit__(self, et, ev, tb): pass + def _log_intro(self, logger, **log_kwargs): + """ + Log PyROS solver introductory messages. + + Parameters + ---------- + logger : logging.Logger + Logger through which to emit messages. + **log_kwargs : dict, optional + Keyword arguments to ``logger.log()`` callable. + Should not include `msg`. + """ + logger.log(msg="=" * self._LOG_LINE_LENGTH, **log_kwargs) + logger.log( + msg=f"PyROS: The Pyomo Robust Optimization Solver, v{self.version()}.", + **log_kwargs, + ) + + # git_info_str = ", ".join( + # f"{field}: {val}" for field, val in _get_pyomo_git_info().items() + # ) + version_info = _get_pyomo_version_info() + version_info_str = ' ' * len("PyROS: ") + ("\n" + ' ' * len("PyROS: ")).join( + f"{key}: {val}" for key, val in version_info.items() + ) + logger.log(msg=version_info_str, **log_kwargs) + logger.log( + msg=( + f"{' ' * len('PyROS:')} " + f"Invoked at UTC {datetime.utcnow().isoformat()}" + ), + **log_kwargs, + ) + logger.log(msg="", **log_kwargs) + logger.log( + msg=("Developed by: Natalie M. Isenberg (1), Jason A. F. Sherman (1),"), + **log_kwargs, + ) + logger.log( + msg=( + f"{' ' * len('Developed by:')} " + "John D. Siirola (2), Chrysanthos E. Gounaris (1)" + ), + **log_kwargs, + ) + logger.log( + msg=( + "(1) Carnegie Mellon University, " "Department of Chemical Engineering" + ), + **log_kwargs, + ) + logger.log( + msg="(2) Sandia National Laboratories, Center for Computing Research", + **log_kwargs, + ) + logger.log(msg="", **log_kwargs) + logger.log( + msg=( + "The developers gratefully acknowledge support " + "from the U.S. Department" + ), + **log_kwargs, + ) + logger.log( + msg=( + "of Energy's " + "Institute for the Design of Advanced Energy Systems (IDAES)." + ), + **log_kwargs, + ) + logger.log(msg="=" * self._LOG_LINE_LENGTH, **log_kwargs) + + def _log_disclaimer(self, logger, **log_kwargs): + """ + Log PyROS solver disclaimer messages. + + Parameters + ---------- + logger : logging.Logger + Logger through which to emit messages. + **log_kwargs : dict, optional + Keyword arguments to ``logger.log()`` callable. + Should not include `msg`. + """ + disclaimer_header = " DISCLAIMER ".center(self._LOG_LINE_LENGTH, "=") + + logger.log(msg=disclaimer_header, **log_kwargs) + logger.log(msg="PyROS is still under development. ", **log_kwargs) + logger.log( + msg=( + "Please provide feedback and/or report any issues by creating " + "a ticket at" + ), + **log_kwargs, + ) + logger.log(msg="https://github.com/Pyomo/pyomo/issues/new/choose", **log_kwargs) + logger.log(msg="=" * self._LOG_LINE_LENGTH, **log_kwargs) + + def _log_config(self, logger, config, exclude_options=None, **log_kwargs): + """ + Log PyROS solver options. + + Parameters + ---------- + logger : logging.Logger + Logger for the solver options. + config : ConfigDict + PyROS solver options. + exclude_options : None or iterable of str, optional + Options (keys of the ConfigDict) to exclude from + logging. If `None` passed, then the names of the + required arguments to ``self.solve()`` are skipped. + **log_kwargs : dict, optional + Keyword arguments to each statement of ``logger.log()``. + """ + # log solver options + if exclude_options is None: + exclude_options = [ + "first_stage_variables", + "second_stage_variables", + "uncertain_params", + "uncertainty_set", + "local_solver", + "global_solver", + ] + + logger.log(msg="Solver options:", **log_kwargs) + for key, val in config.items(): + if key not in exclude_options: + logger.log(msg=f" {key}={val!r}", **log_kwargs) + logger.log(msg="-" * self._LOG_LINE_LENGTH, **log_kwargs) + def solve( self, model, @@ -742,15 +915,26 @@ def solve( model_data = ROSolveResults() model_data.timing = Bunch() - # === Set up logger for logging results - with time_code(model_data.timing, 'total', is_main_timer=True): - config.progress_logger.setLevel(logging.INFO) - - # === PREAMBLE - output_logger(config=config, preamble=True, version=str(self.version())) + # === Start timer, run the algorithm + model_data.timing = TimingData() + with time_code( + timing_data_obj=model_data.timing, + code_block_name="main", + is_main_timer=True, + ): + # output intro and disclaimer + self._log_intro(logger=config.progress_logger, level=logging.INFO) + self._log_disclaimer(logger=config.progress_logger, level=logging.INFO) + self._log_config( + logger=config.progress_logger, + config=config, + exclude_options=None, + level=logging.INFO, + ) - # === DISCLAIMER - output_logger(config=config, disclaimer=True) + # begin preprocessing + config.progress_logger.info("Preprocessing...") + model_data.timing.start_timer("main.preprocessing") # === A block to hold list-type data to make cloning easy util = Block(concrete=True) @@ -783,6 +967,7 @@ def solve( ) assert len(active_objs) == 1 active_obj = active_objs[0] + active_obj_original_sense = active_obj.sense recast_to_min_obj(model_data.working_model, active_obj) # === Determine first and second-stage objectives @@ -831,10 +1016,18 @@ def solve( if "bound_con" in c.name: wm_util.ssv_bounds.append(c) + model_data.timing.stop_timer("main.preprocessing") + preprocessing_time = model_data.timing.get_total_time("main.preprocessing") + config.progress_logger.info( + f"Done preprocessing; required wall time of " + f"{preprocessing_time:.3f}s." + ) + # === Solve and load solution into model pyros_soln, final_iter_separation_solns = ROSolver_iterative_solve( model_data, config ) + IterationLogRecord.log_header_rule(config.progress_logger.info) return_soln = ROSolveResults() if pyros_soln is not None and final_iter_separation_solns is not None: @@ -846,31 +1039,24 @@ def solve( ): load_final_solution(model_data, pyros_soln.master_soln, config) - # === Return time info - model_data.total_cpu_time = get_main_elapsed_time(model_data.timing) - iterations = pyros_soln.total_iters + 1 - - # === Return config to user - return_soln.config = config - # Report the negative of the objective value if it was originally maximize, since we use the minimize form in the algorithm - if next(model.component_data_objects(Objective)).sense == maximize: - negation = -1 - else: - negation = 1 + # account for sense of the original model objective + # when reporting the final PyROS (master) objective, + # since maximization objective is changed to + # minimization objective during preprocessing if config.objective_focus == ObjectiveType.nominal: - return_soln.final_objective_value = negation * value( - pyros_soln.master_soln.master_model.obj + return_soln.final_objective_value = ( + active_obj_original_sense + * value(pyros_soln.master_soln.master_model.obj) ) elif config.objective_focus == ObjectiveType.worst_case: - return_soln.final_objective_value = negation * value( - pyros_soln.master_soln.master_model.zeta + return_soln.final_objective_value = ( + active_obj_original_sense + * value(pyros_soln.master_soln.master_model.zeta) ) return_soln.pyros_termination_condition = ( pyros_soln.pyros_termination_condition ) - - return_soln.time = model_data.total_cpu_time - return_soln.iterations = iterations + return_soln.iterations = pyros_soln.total_iters + 1 # === Remove util block model.del_component(model_data.util_block) @@ -878,12 +1064,25 @@ def solve( del pyros_soln.util_block del pyros_soln.working_model else: + return_soln.final_objective_value = None return_soln.pyros_termination_condition = ( pyrosTerminationCondition.robust_infeasible ) - return_soln.final_objective_value = None - return_soln.time = get_main_elapsed_time(model_data.timing) return_soln.iterations = 0 + + return_soln.config = config + return_soln.time = model_data.timing.get_total_time("main") + + # log termination-related messages + config.progress_logger.info(return_soln.pyros_termination_condition.message) + config.progress_logger.info("-" * self._LOG_LINE_LENGTH) + config.progress_logger.info(f"Timing breakdown:\n\n{model_data.timing}") + config.progress_logger.info("-" * self._LOG_LINE_LENGTH) + config.progress_logger.info(return_soln) + config.progress_logger.info("-" * self._LOG_LINE_LENGTH) + config.progress_logger.info("All done. Exiting PyROS.") + config.progress_logger.info("=" * self._LOG_LINE_LENGTH) + return return_soln diff --git a/pyomo/contrib/pyros/pyros_algorithm_methods.py b/pyomo/contrib/pyros/pyros_algorithm_methods.py index 7a0c990d549..df0d539a70d 100644 --- a/pyomo/contrib/pyros/pyros_algorithm_methods.py +++ b/pyomo/contrib/pyros/pyros_algorithm_methods.py @@ -11,14 +11,14 @@ ObjectiveType, get_time_from_solver, pyrosTerminationCondition, + IterationLogRecord, ) -from pyomo.contrib.pyros.util import ( - get_main_elapsed_time, - output_logger, - coefficient_matching, -) +from pyomo.contrib.pyros.util import get_main_elapsed_time, coefficient_matching from pyomo.core.base import value -from pyomo.common.collections import ComponentSet +from pyomo.common.collections import ComponentSet, ComponentMap +from pyomo.core.base.var import _VarData as VarData +from itertools import chain +from pyomo.common.dependencies import numpy as np def update_grcs_solve_data( @@ -47,6 +47,150 @@ def update_grcs_solve_data( return +def get_dr_var_to_scaled_expr_map( + decision_rule_eqns, second_stage_vars, uncertain_params, decision_rule_vars +): + """ + Generate mapping from decision rule variables + to their terms in a model's DR expression. + """ + var_to_scaled_expr_map = ComponentMap() + ssv_dr_eq_zip = zip(second_stage_vars, decision_rule_eqns) + for ssv_idx, (ssv, dr_eq) in enumerate(ssv_dr_eq_zip): + for term in dr_eq.body.args: + is_ssv_term = ( + isinstance(term.args[0], int) + and term.args[0] == -1 + and isinstance(term.args[1], VarData) + ) + if not is_ssv_term: + dr_var = term.args[1] + var_to_scaled_expr_map[dr_var] = term + + return var_to_scaled_expr_map + + +def evaluate_and_log_component_stats(model_data, separation_model, config): + """ + Evaluate and log model component statistics. + """ + IterationLogRecord.log_header_rule(config.progress_logger.info) + config.progress_logger.info("Model statistics:") + # print model statistics + dr_var_set = ComponentSet( + chain( + *tuple( + indexed_dr_var.values() + for indexed_dr_var in model_data.working_model.util.decision_rule_vars + ) + ) + ) + first_stage_vars = [ + var + for var in model_data.working_model.util.first_stage_variables + if var not in dr_var_set + ] + + # account for epigraph constraint + sep_model_epigraph_con = getattr(separation_model, "epigraph_constr", None) + has_epigraph_con = sep_model_epigraph_con is not None + + num_fsv = len(first_stage_vars) + num_ssv = len(model_data.working_model.util.second_stage_variables) + num_sv = len(model_data.working_model.util.state_vars) + num_dr_vars = len(dr_var_set) + num_vars = int(has_epigraph_con) + num_fsv + num_ssv + num_sv + num_dr_vars + + num_uncertain_params = len(model_data.working_model.util.uncertain_params) + + eq_cons = [ + con + for con in model_data.working_model.component_data_objects( + Constraint, active=True + ) + if con.equality + ] + dr_eq_set = ComponentSet( + chain( + *tuple( + indexed_dr_eq.values() + for indexed_dr_eq in model_data.working_model.util.decision_rule_eqns + ) + ) + ) + num_eq_cons = len(eq_cons) + num_dr_cons = len(dr_eq_set) + num_coefficient_matching_cons = len( + getattr(model_data.working_model, "coefficient_matching_constraints", []) + ) + num_other_eq_cons = num_eq_cons - num_dr_cons - num_coefficient_matching_cons + + # get performance constraints as referenced in the separation + # model object + new_sep_con_map = separation_model.util.map_new_constraint_list_to_original_con + perf_con_set = ComponentSet( + new_sep_con_map.get(con, con) + for con in separation_model.util.performance_constraints + ) + is_epigraph_con_first_stage = ( + has_epigraph_con and sep_model_epigraph_con not in perf_con_set + ) + working_model_perf_con_set = ComponentSet( + model_data.working_model.find_component(new_sep_con_map.get(con, con)) + for con in separation_model.util.performance_constraints + if con is not None + ) + + num_perf_cons = len(separation_model.util.performance_constraints) + num_fsv_bounds = sum( + int(var.lower is not None) + int(var.upper is not None) + for var in first_stage_vars + ) + ineq_con_set = [ + con + for con in model_data.working_model.component_data_objects( + Constraint, active=True + ) + if not con.equality + ] + num_fsv_ineqs = ( + num_fsv_bounds + + len([con for con in ineq_con_set if con not in working_model_perf_con_set]) + + is_epigraph_con_first_stage + ) + num_ineq_cons = len(ineq_con_set) + has_epigraph_con + num_fsv_bounds + + config.progress_logger.info(f"{' Number of variables'} : {num_vars}") + config.progress_logger.info(f"{' Epigraph variable'} : {int(has_epigraph_con)}") + config.progress_logger.info(f"{' First-stage variables'} : {num_fsv}") + config.progress_logger.info(f"{' Second-stage variables'} : {num_ssv}") + config.progress_logger.info(f"{' State variables'} : {num_sv}") + config.progress_logger.info(f"{' Decision rule variables'} : {num_dr_vars}") + config.progress_logger.info( + f"{' Number of uncertain parameters'} : {num_uncertain_params}" + ) + config.progress_logger.info( + f"{' Number of constraints'} : " f"{num_ineq_cons + num_eq_cons}" + ) + config.progress_logger.info(f"{' Equality constraints'} : {num_eq_cons}") + config.progress_logger.info( + f"{' Coefficient matching constraints'} : " + f"{num_coefficient_matching_cons}" + ) + config.progress_logger.info(f"{' Decision rule equations'} : {num_dr_cons}") + config.progress_logger.info( + f"{' All other equality constraints'} : " f"{num_other_eq_cons}" + ) + config.progress_logger.info(f"{' Inequality constraints'} : {num_ineq_cons}") + config.progress_logger.info( + f"{' First-stage inequalities (incl. certain var bounds)'} : " + f"{num_fsv_ineqs}" + ) + config.progress_logger.info( + f"{' Performance constraints (incl. var bounds)'} : {num_perf_cons}" + ) + + def ROSolver_iterative_solve(model_data, config): ''' GRCS algorithm implementation @@ -75,20 +219,23 @@ def ROSolver_iterative_solve(model_data, config): config=config, ) if not coeff_matching_success and not robust_infeasible: - raise ValueError( - "Equality constraint \"%s\" cannot be guaranteed to be robustly feasible, " - "given the current partitioning between first-stage, second-stage and state variables. " - "You might consider editing this constraint to reference some second-stage " - "and/or state variable(s)." % c.name + config.progress_logger.error( + f"Equality constraint {c.name!r} cannot be guaranteed to " + "be robustly feasible, given the current partitioning " + "among first-stage, second-stage, and state variables. " + "Consider editing this constraint to reference some " + "second-stage and/or state variable(s)." ) + raise ValueError("Coefficient matching unsuccessful. See the solver logs.") elif not coeff_matching_success and robust_infeasible: config.progress_logger.info( "PyROS has determined that the model is robust infeasible. " - "One reason for this is that equality constraint \"%s\" cannot be satisfied " - "against all realizations of uncertainty, " - "given the current partitioning between first-stage, second-stage and state variables. " - "You might consider editing this constraint to reference some (additional) second-stage " - "and/or state variable(s)." % c.name + f"One reason for this is that the equality constraint {c.name} " + "cannot be satisfied against all realizations of uncertainty, " + "given the current partitioning between " + "first-stage, second-stage, and state variables. " + "Consider editing this constraint to reference some (additional) " + "second-stage and/or state variable(s)." ) return None, None else: @@ -156,6 +303,10 @@ def ROSolver_iterative_solve(model_data, config): model_data=master_data, config=config ) + evaluate_and_log_component_stats( + model_data=model_data, separation_model=separation_model, config=config + ) + # === Create separation problem data container object and add information to catalog during solve separation_data = SeparationProblemData() separation_data.separation_model = separation_model @@ -204,6 +355,34 @@ def ROSolver_iterative_solve(model_data, config): dr_var_lists_original = [] dr_var_lists_polished = [] + # set up first-stage variable and DR variable sets + master_dr_var_set = ComponentSet( + chain( + *tuple( + indexed_var.values() + for indexed_var in master_data.master_model.scenarios[ + 0, 0 + ].util.decision_rule_vars + ) + ) + ) + master_fsv_set = ComponentSet( + var + for var in master_data.master_model.scenarios[0, 0].util.first_stage_variables + if var not in master_dr_var_set + ) + previous_master_fsv_vals = ComponentMap((var, None) for var in master_fsv_set) + previous_master_dr_var_vals = ComponentMap((var, None) for var in master_dr_var_set) + + nom_master_util_blk = master_data.master_model.scenarios[0, 0].util + dr_var_scaled_expr_map = get_dr_var_to_scaled_expr_map( + decision_rule_vars=nom_master_util_blk.decision_rule_vars, + decision_rule_eqns=nom_master_util_blk.decision_rule_eqns, + second_stage_vars=nom_master_util_blk.second_stage_variables, + uncertain_params=nom_master_util_blk.uncertain_params, + ) + + IterationLogRecord.log_header(config.progress_logger.info) k = 0 master_statuses = [] while config.max_iter == -1 or k < config.max_iter: @@ -216,7 +395,7 @@ def ROSolver_iterative_solve(model_data, config): ) # === Solve Master Problem - config.progress_logger.info("PyROS working on iteration %s..." % k) + config.progress_logger.debug(f"PyROS working on iteration {k}...") master_soln = master_problem_methods.solve_master( model_data=master_data, config=config ) @@ -239,7 +418,6 @@ def ROSolver_iterative_solve(model_data, config): is pyrosTerminationCondition.robust_infeasible ): term_cond = pyrosTerminationCondition.robust_infeasible - output_logger(config=config, robust_infeasible=True) elif ( master_soln.pyros_termination_condition is pyrosTerminationCondition.subsolver_error @@ -257,6 +435,19 @@ def ROSolver_iterative_solve(model_data, config): pyrosTerminationCondition.time_out, pyrosTerminationCondition.robust_infeasible, }: + log_record = IterationLogRecord( + iteration=k, + objective=None, + first_stage_var_shift=None, + dr_var_shift=None, + num_violated_cons=None, + max_violation=None, + dr_polishing_success=None, + all_sep_problems_solved=None, + global_separation=None, + elapsed_time=get_main_elapsed_time(model_data.timing), + ) + log_record.log(config.progress_logger.info) update_grcs_solve_data( pyros_soln=model_data, k=k, @@ -280,6 +471,7 @@ def ROSolver_iterative_solve(model_data, config): nominal_data.nom_second_stage_cost = master_soln.second_stage_objective nominal_data.nom_obj = value(master_data.master_model.obj) + polishing_successful = True if ( config.decision_rule_order != 0 and len(config.second_stage_variables) > 0 @@ -294,7 +486,10 @@ def ROSolver_iterative_solve(model_data, config): vals.append(dvar.value) dr_var_lists_original.append(vals) - polishing_results = master_problem_methods.minimize_dr_vars( + ( + polishing_results, + polishing_successful, + ) = master_problem_methods.minimize_dr_vars( model_data=master_data, config=config ) timing_data.total_dr_polish_time += get_time_from_solver(polishing_results) @@ -308,11 +503,55 @@ def ROSolver_iterative_solve(model_data, config): vals.append(dvar.value) dr_var_lists_polished.append(vals) - # === Check if time limit reached - elapsed = get_main_elapsed_time(model_data.timing) + # get current first-stage and DR variable values + # and compare with previous first-stage and DR variable + # values + current_master_fsv_vals = ComponentMap( + (var, value(var)) for var in master_fsv_set + ) + current_master_dr_var_vals = ComponentMap( + (var, value(var)) for var, expr in dr_var_scaled_expr_map.items() + ) + if k > 0: + first_stage_var_shift = ( + max( + abs(current_master_fsv_vals[var] - previous_master_fsv_vals[var]) + for var in previous_master_fsv_vals + ) + if current_master_fsv_vals + else None + ) + dr_var_shift = ( + max( + abs( + current_master_dr_var_vals[var] + - previous_master_dr_var_vals[var] + ) + for var in previous_master_dr_var_vals + ) + if current_master_dr_var_vals + else None + ) + else: + first_stage_var_shift = None + dr_var_shift = None + + # === Check if time limit reached after polishing if config.time_limit: + elapsed = get_main_elapsed_time(model_data.timing) if elapsed >= config.time_limit: - output_logger(config=config, time_out=True, elapsed=elapsed) + iter_log_record = IterationLogRecord( + iteration=k, + objective=value(master_data.master_model.obj), + first_stage_var_shift=first_stage_var_shift, + dr_var_shift=dr_var_shift, + num_violated_cons=None, + max_violation=None, + dr_polishing_success=polishing_successful, + all_sep_problems_solved=None, + global_separation=None, + elapsed_time=elapsed, + ) update_grcs_solve_data( pyros_soln=model_data, k=k, @@ -322,6 +561,7 @@ def ROSolver_iterative_solve(model_data, config): separation_data=separation_data, master_soln=master_soln, ) + iter_log_record.log(config.progress_logger.info) return model_data, [] # === Set up for the separation problem @@ -376,10 +616,39 @@ def ROSolver_iterative_solve(model_data, config): separation_results.violating_param_realization ) + scaled_violations = [ + solve_call_res.scaled_violations[con] + for con, solve_call_res in separation_results.main_loop_results.solver_call_results.items() + if solve_call_res.scaled_violations is not None + ] + if scaled_violations: + max_sep_con_violation = max(scaled_violations) + else: + max_sep_con_violation = None + num_violated_cons = len(separation_results.violated_performance_constraints) + + all_sep_problems_solved = ( + len(scaled_violations) == len(separation_model.util.performance_constraints) + and not separation_results.subsolver_error + and not separation_results.time_out + ) + + iter_log_record = IterationLogRecord( + iteration=k, + objective=value(master_data.master_model.obj), + first_stage_var_shift=first_stage_var_shift, + dr_var_shift=dr_var_shift, + num_violated_cons=num_violated_cons, + max_violation=max_sep_con_violation, + dr_polishing_success=polishing_successful, + all_sep_problems_solved=all_sep_problems_solved, + global_separation=separation_results.solved_globally, + elapsed_time=get_main_elapsed_time(model_data.timing), + ) + # terminate on time limit elapsed = get_main_elapsed_time(model_data.timing) if separation_results.time_out: - output_logger(config=config, time_out=True, elapsed=elapsed) termination_condition = pyrosTerminationCondition.time_out update_grcs_solve_data( pyros_soln=model_data, @@ -390,6 +659,7 @@ def ROSolver_iterative_solve(model_data, config): separation_data=separation_data, master_soln=master_soln, ) + iter_log_record.log(config.progress_logger.info) return model_data, separation_results # terminate on separation subsolver error @@ -404,24 +674,26 @@ def ROSolver_iterative_solve(model_data, config): separation_data=separation_data, master_soln=master_soln, ) + iter_log_record.log(config.progress_logger.info) return model_data, separation_results # === Check if we terminate due to robust optimality or feasibility, # or in the event of bypassing global separation, no violations robustness_certified = separation_results.robustness_certified if robustness_certified: - output_logger( - config=config, bypass_global_separation=config.bypass_global_separation - ) + if config.bypass_global_separation: + config.progress_logger.warning( + "Option to bypass global separation was chosen. " + "Robust feasibility and optimality of the reported " + "solution are not guaranteed." + ) robust_optimal = ( config.solve_master_globally and config.objective_focus is ObjectiveType.worst_case ) if robust_optimal: - output_logger(config=config, robust_optimal=True) termination_condition = pyrosTerminationCondition.robust_optimal else: - output_logger(config=config, robust_feasible=True) termination_condition = pyrosTerminationCondition.robust_feasible update_grcs_solve_data( pyros_soln=model_data, @@ -432,6 +704,7 @@ def ROSolver_iterative_solve(model_data, config): separation_data=separation_data, master_soln=master_soln, ) + iter_log_record.log(config.progress_logger.info) return model_data, separation_results # === Add block to master at violation @@ -443,13 +716,21 @@ def ROSolver_iterative_solve(model_data, config): separation_results.violating_param_realization ) + config.progress_logger.debug("Points added to master:") + config.progress_logger.debug( + np.array([pt for pt in separation_data.points_added_to_master]) + ) + k += 1 + iter_log_record.log(config.progress_logger.info) + previous_master_fsv_vals = current_master_fsv_vals + previous_master_dr_var_vals = current_master_dr_var_vals + # Iteration limit reached - output_logger(config=config, max_iter=True) update_grcs_solve_data( pyros_soln=model_data, - k=k, + k=k - 1, # remove last increment to fix iteration count term_cond=pyrosTerminationCondition.max_iter, nominal_data=nominal_data, timing_data=timing_data, diff --git a/pyomo/contrib/pyros/separation_problem_methods.py b/pyomo/contrib/pyros/separation_problem_methods.py index 2c41c869474..61f347b418d 100644 --- a/pyomo/contrib/pyros/separation_problem_methods.py +++ b/pyomo/contrib/pyros/separation_problem_methods.py @@ -6,7 +6,7 @@ from pyomo.core.base import Var, Param from pyomo.common.collections import ComponentSet, ComponentMap from pyomo.common.dependencies import numpy as np -from pyomo.contrib.pyros.util import ObjectiveType, get_time_from_solver, output_logger +from pyomo.contrib.pyros.util import ObjectiveType, get_time_from_solver from pyomo.contrib.pyros.solve_data import ( DiscreteSeparationSolveCallResults, SeparationSolveCallResults, @@ -536,6 +536,54 @@ def get_worst_discrete_separation_solution( ) +def get_con_name_repr(separation_model, con, with_orig_name=True, with_obj_name=True): + """ + Get string representation of performance constraint + and any other modeling components to which it has + been mapped. + + Parameters + ---------- + separation_model : ConcreteModel + Separation model. + con : ScalarConstraint or ConstraintData + Constraint for which to get the representation. + with_orig_name : bool, optional + If constraint was added during construction of the + separation problem (i.e. if the constraint is a member of + in `separation_model.util.new_constraints`), + include the name of the original constraint from which + `perf_con` was created. + with_obj_name : bool, optional + Include name of separation model objective to which + constraint is mapped. Applicable only to performance + constraints of the separation problem. + + Returns + ------- + str + Constraint name representation. + """ + + qual_strs = [] + if with_orig_name: + # check performance constraint was not added + # at construction of separation problem + orig_con = separation_model.util.map_new_constraint_list_to_original_con.get( + con, con + ) + if orig_con is not con: + qual_strs.append(f"originally {orig_con.name!r}") + if with_obj_name: + objectives_map = separation_model.util.map_obj_to_constr + separation_obj = objectives_map[con] + qual_strs.append(f"mapped to objective {separation_obj.name!r}") + + final_qual_str = f" ({', '.join(qual_strs)})" if qual_strs else "" + + return f"{con.name!r}{final_qual_str}" + + def perform_separation_loop(model_data, config, solve_globally): """ Loop through, and solve, PyROS separation problems to @@ -633,12 +681,22 @@ def perform_separation_loop(model_data, config, solve_globally): ) all_solve_call_results = ComponentMap() - for priority, perf_constraints in sorted_priority_groups.items(): + priority_groups_enum = enumerate(sorted_priority_groups.items()) + for group_idx, (priority, perf_constraints) in priority_groups_enum: priority_group_solve_call_results = ComponentMap() - for perf_con in perf_constraints: - # config.progress_logger.info( - # f"Separating constraint {perf_con.name}" - # ) + for idx, perf_con in enumerate(perf_constraints): + # log progress of separation loop + solve_adverb = "Globally" if solve_globally else "Locally" + config.progress_logger.debug( + f"{solve_adverb} separating performance constraint " + f"{get_con_name_repr(model_data.separation_model, perf_con)} " + f"(priority {priority}, priority group {group_idx + 1} of " + f"{len(sorted_priority_groups)}, " + f"constraint {idx + 1} of {len(perf_constraints)} " + "in priority group, " + f"{len(all_solve_call_results) + idx + 1} of " + f"{len(all_performance_constraints)} total)" + ) # solve separation problem for this performance constraint if uncertainty_set_is_discrete: @@ -689,21 +747,33 @@ def perform_separation_loop(model_data, config, solve_globally): ) # # auxiliary log messages - # objectives_map = ( - # model_data.separation_model.util.map_obj_to_constr - # ) - # violated_con_name = list(objectives_map.keys())[ - # worst_case_perf_con - # ] - # config.progress_logger.info( - # f"Violation found for constraint {violated_con_name} " - # "under realization " - # f"{worst_case_res.violating_param_realization}" - # ) + violated_con_names = "\n ".join( + get_con_name_repr(model_data.separation_model, con) + for con, res in all_solve_call_results.items() + if res.found_violation + ) + config.progress_logger.debug( + f"Violated constraints:\n {violated_con_names} " + ) + config.progress_logger.debug( + "Worst-case constraint: " + f"{get_con_name_repr(model_data.separation_model, worst_case_perf_con)} " + "under realization " + f"{worst_case_res.violating_param_realization}." + ) + config.progress_logger.debug( + f"Maximal scaled violation " + f"{worst_case_res.scaled_violations[worst_case_perf_con]} " + "from this constraint " + "exceeds the robust feasibility tolerance " + f"{config.robust_feasibility_tolerance}" + ) # violating separation problem solution now chosen. # exit loop break + else: + config.progress_logger.debug("No violated performance constraints found.") return SeparationLoopResults( solver_call_results=all_solve_call_results, @@ -786,7 +856,7 @@ def evaluate_performance_constraint_violations( return (violating_param_realization, scaled_violations, constraint_violated) -def initialize_separation(model_data, config): +def initialize_separation(perf_con_to_maximize, model_data, config): """ Initialize separation problem variables, and fix all first-stage variables to their corresponding values from most recent @@ -794,6 +864,9 @@ def initialize_separation(model_data, config): Parameters ---------- + perf_con_to_maximize : ConstraintData + Performance constraint whose violation is to be maximized + for the separation problem of interest. model_data : SeparationProblemData Separation problem data. config : ConfigDict @@ -875,13 +948,34 @@ def get_parent_master_blk(var): "All h(x,q) type constraints must be deactivated in separation." ) - # check: initial point feasible? + # confirm the initial point is feasible for cases where + # we expect it to be (i.e. non-discrete uncertainty sets). + # otherwise, log the violated constraints + tol = ABS_CON_CHECK_FEAS_TOL + perf_con_name_repr = get_con_name_repr( + separation_model=model_data.separation_model, + con=perf_con_to_maximize, + with_orig_name=True, + with_obj_name=True, + ) + uncertainty_set_is_discrete = ( + config.uncertainty_set.geometry is Geometry.DISCRETE_SCENARIOS + ) for con in sep_model.component_data_objects(Constraint, active=True): - lb, val, ub = value(con.lb), value(con.body), value(con.ub) - lb_viol = val < lb - ABS_CON_CHECK_FEAS_TOL if lb is not None else False - ub_viol = val > ub + ABS_CON_CHECK_FEAS_TOL if ub is not None else False - if lb_viol or ub_viol: - config.progress_logger.debug(con.name, lb, val, ub) + lslack, uslack = con.lslack(), con.uslack() + if (lslack < -tol or uslack < -tol) and not uncertainty_set_is_discrete: + con_name_repr = get_con_name_repr( + separation_model=model_data.separation_model, + con=con, + with_orig_name=True, + with_obj_name=False, + ) + config.progress_logger.debug( + f"Initial point for separation of performance constraint " + f"{perf_con_name_repr} violates the model constraint " + f"{con_name_repr} by more than {tol}. " + f"(lslack={con.lslack()}, uslack={con.uslack()})" + ) locally_acceptable = {tc.optimal, tc.locallyOptimal, tc.globallyOptimal} @@ -929,8 +1023,17 @@ def solver_call_separation( solver_status_dict = {} nlp_model = model_data.separation_model + # get name of constraint for loggers + con_name_repr = get_con_name_repr( + separation_model=nlp_model, + con=perf_con_to_maximize, + with_orig_name=True, + with_obj_name=True, + ) + solve_mode = "global" if solve_globally else "local" + # === Initialize separation problem; fix first-stage variables - initialize_separation(model_data, config) + initialize_separation(perf_con_to_maximize, model_data, config) separation_obj.activate() @@ -942,10 +1045,18 @@ def solver_call_separation( subsolver_error=False, ) timer = TicTocTimer() - for opt in solvers: + for idx, opt in enumerate(solvers): + if idx > 0: + config.progress_logger.warning( + f"Invoking backup solver {opt!r} " + f"(solver {idx + 1} of {len(solvers)}) for {solve_mode} " + f"separation of performance constraint {con_name_repr} " + f"in iteration {model_data.iteration}." + ) orig_setting, custom_setting_present = adjust_solver_time_settings( model_data.timing, opt, config ) + model_data.timing.start_timer(f"main.{solve_mode}_separation") timer.tic(msg=None) try: results = opt.solve( @@ -958,14 +1069,17 @@ def solver_call_separation( # account for possible external subsolver errors # (such as segmentation faults, function evaluation # errors, etc.) + adverb = "globally" if solve_globally else "locally" config.progress_logger.error( - f"Solver {repr(opt)} encountered exception attempting to " - "optimize separation problem in iteration " - f"{model_data.iteration}" + f"Optimizer {repr(opt)} ({idx + 1} of {len(solvers)}) " + f"encountered exception attempting " + f"to {adverb} solve separation problem for constraint " + f"{con_name_repr} in iteration {model_data.iteration}." ) raise else: setattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR, timer.toc(msg=None)) + model_data.timing.stop_timer(f"main.{solve_mode}_separation") finally: revert_solver_max_time_adjustment( opt, orig_setting, custom_setting_present, config @@ -1014,15 +1128,25 @@ def solver_call_separation( separation_obj.deactivate() return solve_call_results + else: + config.progress_logger.debug( + f"Solver {opt} ({idx + 1} of {len(solvers)}) " + f"failed for {solve_mode} separation of performance " + f"constraint {con_name_repr} in iteration " + f"{model_data.iteration}. Termination condition: " + f"{results.solver.termination_condition!r}." + ) + config.progress_logger.debug(f"Results:\n{results.solver}") # All subordinate solvers failed to optimize model to appropriate # termination condition. PyROS will terminate with subsolver # error. At this point, export model if desired solve_call_results.subsolver_error = True save_dir = config.subproblem_file_directory + serialization_msg = "" if save_dir and config.keepfiles: objective = separation_obj.name - name = os.path.join( + output_problem_path = os.path.join( save_dir, ( config.uncertainty_set.type @@ -1035,15 +1159,23 @@ def solver_call_separation( + ".bar" ), ) - nlp_model.write(name, io_options={'symbolic_solver_labels': True}) - output_logger( - config=config, - separation_error=True, - filename=name, - iteration=model_data.iteration, - objective=objective, - status_dict=solver_status_dict, + nlp_model.write( + output_problem_path, io_options={'symbolic_solver_labels': True} ) + serialization_msg = ( + " For debugging, problem has been serialized to the file " + f"{output_problem_path!r}." + ) + solve_call_results.message = ( + "Could not successfully solve separation problem of iteration " + f"{model_data.iteration} " + f"for performance constraint {con_name_repr} with any of the " + f"provided subordinate {solve_mode} optimizers. " + f"(Termination statuses: " + f"{[str(term_cond) for term_cond in solver_status_dict.values()]}.)" + f"{serialization_msg}" + ) + config.progress_logger.warning(solve_call_results.message) separation_obj.deactivate() diff --git a/pyomo/contrib/pyros/solve_data.py b/pyomo/contrib/pyros/solve_data.py index 63e7fdd7ebd..40a52757bae 100644 --- a/pyomo/contrib/pyros/solve_data.py +++ b/pyomo/contrib/pyros/solve_data.py @@ -5,17 +5,71 @@ class ROSolveResults(object): """ - Container for solve-instance data returned to the user after solving with PyROS. + PyROS solver results object. - Attributes: - :pyros_termination_condition: termination condition of the PyROS algorithm - :config: the config block for this solve instance - :time: Total solver CPU time - :iterations: total iterations done by PyROS solver - :final_objective_value: objective function value at termination + Parameters + ---------- + config : ConfigDict, optional + User-specified solver settings. + iterations : int, optional + Number of iterations required. + time : float, optional + Total elapsed time (or wall time), in seconds. + final_objective_value : float, optional + Final objective function value to report. + pyros_termination_condition : pyrosTerminationCondition, optional + PyROS-specific termination condition. + + Attributes + ---------- + config : ConfigDict, optional + User-specified solver settings. + iterations : int, optional + Number of iterations required by PyROS. + time : float, optional + Total elapsed time (or wall time), in seconds. + final_objective_value : float, optional + Final objective function value to report. + pyros_termination_condition : pyros.util.pyrosTerminationStatus + Indicator of the manner of termination. """ - pass + def __init__( + self, + config=None, + iterations=None, + time=None, + final_objective_value=None, + pyros_termination_condition=None, + ): + """Initialize self (see class docstring).""" + self.config = config + self.iterations = iterations + self.time = time + self.final_objective_value = final_objective_value + self.pyros_termination_condition = pyros_termination_condition + + def __str__(self): + """ + Generate string representation of self. + Does not include any information about `self.config`. + """ + lines = ["Termination stats:"] + attr_name_format_dict = { + "iterations": ("Iterations", "f'{val}'"), + "time": ("Solve time (wall s)", "f'{val:.3f}'"), + "final_objective_value": ("Final objective value", "f'{val:.4e}'"), + "pyros_termination_condition": ("Termination condition", "f'{val}'"), + } + attr_desc_pad_length = max( + len(desc) for desc, _ in attr_name_format_dict.values() + ) + for attr_name, (attr_desc, fmt_str) in attr_name_format_dict.items(): + val = getattr(self, attr_name) + val_str = eval(fmt_str) if val is not None else str(val) + lines.append(f" {attr_desc:<{attr_desc_pad_length}s} : {val_str}") + + return "\n".join(lines) class MasterProblemData(object): @@ -443,6 +497,7 @@ class SeparationResults: ---------- local_separation_loop_results global_separation_loop_results + main_loop_results subsolver_error time_out solved_locally @@ -462,7 +517,7 @@ def __init__(self, local_separation_loop_results, global_separation_loop_results @property def time_out(self): """ - Return True if time out found for local or global + bool : True if time out found for local or global separation loop, False otherwise. """ local_time_out = ( @@ -476,7 +531,7 @@ def time_out(self): @property def subsolver_error(self): """ - Return True if subsolver error found for local or global + bool : True if subsolver error found for local or global separation loop, False otherwise. """ local_subsolver_error = ( @@ -490,7 +545,7 @@ def subsolver_error(self): @property def solved_locally(self): """ - Return true if local separation loop was invoked, + bool : true if local separation loop was invoked, False otherwise. """ return self.local_separation_loop_results is not None @@ -498,13 +553,18 @@ def solved_locally(self): @property def solved_globally(self): """ - Return True if global separation loop was invoked, + bool : True if global separation loop was invoked, False otherwise. """ return self.global_separation_loop_results is not None def get_violating_attr(self, attr_name): """ + If separation problems solved globally, returns + value of attribute of global separation loop results. + + Otherwise, if separation problems solved locally, + returns value of attribute of local separation loop results. If local separation loop results specified, return value of attribute of local separation loop results. @@ -526,27 +586,33 @@ def get_violating_attr(self, attr_name): object Attribute value. """ - if self.solved_locally: - local_loop_val = getattr(self.local_separation_loop_results, attr_name) - else: - local_loop_val = None + return getattr(self.main_loop_results, attr_name, None) - if local_loop_val is not None: - attr_val = local_loop_val - else: - if self.solved_globally: - attr_val = getattr(self.global_separation_loop_results, attr_name) - else: - attr_val = None + @property + def worst_case_perf_con(self): + """ + ConstraintData : Performance constraint corresponding to the + separation solution chosen for the next master problem. + """ + return self.get_violating_attr("worst_case_perf_con") - return attr_val + @property + def main_loop_results(self): + """ + SeparationLoopResults : Main separation loop results. + In particular, this is considered to be the global + loop result if solved globally, and the local loop + results otherwise. + """ + if self.solved_globally: + return self.global_separation_loop_results + return self.local_separation_loop_results @property def found_violation(self): """ - bool: True if ``found_violation`` attribute for - local or global separation loop results found - to be True, False otherwise. + bool : True if ``found_violation`` attribute for + main separation loop results is True, False otherwise. """ found_viol = self.get_violating_attr("found_violation") if found_viol is None: diff --git a/pyomo/contrib/pyros/tests/test_grcs.py b/pyomo/contrib/pyros/tests/test_grcs.py index 0d24b799b99..2be73826f61 100644 --- a/pyomo/contrib/pyros/tests/test_grcs.py +++ b/pyomo/contrib/pyros/tests/test_grcs.py @@ -19,6 +19,8 @@ ObjectiveType, pyrosTerminationCondition, coefficient_matching, + TimingData, + IterationLogRecord, ) from pyomo.contrib.pyros.util import replace_uncertain_bounds_with_constraints from pyomo.contrib.pyros.util import get_vars_from_component @@ -33,7 +35,7 @@ solve_master, minimize_dr_vars, ) -from pyomo.contrib.pyros.solve_data import MasterProblemData +from pyomo.contrib.pyros.solve_data import MasterProblemData, ROSolveResults from pyomo.common.dependencies import numpy as np, numpy_available from pyomo.common.dependencies import scipy as sp, scipy_available from pyomo.environ import maximize as pyo_max @@ -59,6 +61,10 @@ sqrt, value, ) +import logging + + +logger = logging.getLogger(__name__) if not (numpy_available and scipy_available): @@ -3567,7 +3573,7 @@ def test_solve_master(self): expr=master_data.master_model.scenarios[0, 0].x ) master_data.iteration = 0 - master_data.timing = Bunch() + master_data.timing = TimingData() box_set = BoxSet(bounds=[(0, 2)]) solver = SolverFactory(global_solver) @@ -3589,8 +3595,11 @@ def test_solve_master(self): ) config.declare("subproblem_file_directory", ConfigValue(default=None)) config.declare("time_limit", ConfigValue(default=None)) + config.declare( + "progress_logger", ConfigValue(default=logging.getLogger(__name__)) + ) - with time_code(master_data.timing, "total", is_main_timer=True): + with time_code(master_data.timing, "main", is_main_timer=True): master_soln = solve_master(master_data, config) self.assertEqual( master_soln.termination_condition, @@ -3866,7 +3875,7 @@ def test_minimize_dr_norm(self): m.working_model.util.state_vars = [] m.working_model.util.first_stage_variables = [] - config = Block() + config = Bunch() config.decision_rule_order = 1 config.objective_focus = ObjectiveType.nominal config.global_solver = SolverFactory('baron') @@ -3874,6 +3883,7 @@ def test_minimize_dr_norm(self): config.tee = False config.solve_master_globally = True config.time_limit = None + config.progress_logger = logging.getLogger(__name__) add_decision_rule_variables(model_data=m, config=config) add_decision_rule_constraints(model_data=m, config=config) @@ -3892,15 +3902,19 @@ def test_minimize_dr_norm(self): master_data.master_model = master master_data.master_model.const_efficiency_applied = False master_data.master_model.linear_efficiency_applied = False + master_data.iteration = 0 - master_data.timing = Bunch() - with time_code(master_data.timing, "total", is_main_timer=True): - results = minimize_dr_vars(model_data=master_data, config=config) + master_data.timing = TimingData() + with time_code(master_data.timing, "main", is_main_timer=True): + results, success = minimize_dr_vars(model_data=master_data, config=config) self.assertEqual( results.solver.termination_condition, TerminationCondition.optimal, msg="Minimize dr norm did not solve to optimality.", ) + self.assertTrue( + success, msg=f"DR polishing success {success}, expected True." + ) @unittest.skipUnless( baron_license_is_valid, "Global NLP solver is not available and licensed." @@ -3957,7 +3971,8 @@ def test_identifying_violating_param_realization(self): baron_license_is_valid, "Global NLP solver is not available and licensed." ) @unittest.skipUnless( - baron_version < (23, 1, 5), "Test known to fail beginning with Baron 23.1.5" + baron_version < (23, 1, 5) or baron_version >= (23, 6, 23), + "Test known to fail for BARON 23.1.5 and versions preceding 23.6.23", ) def test_terminate_with_max_iter(self): m = ConcreteModel() @@ -4004,6 +4019,15 @@ def test_terminate_with_max_iter(self): msg="Returned termination condition is not return max_iter.", ) + self.assertEqual( + results.iterations, + 1, + msg=( + f"Number of iterations in results object is {results.iterations}, " + f"but expected value 1." + ), + ) + @unittest.skipUnless( baron_license_is_valid, "Global NLP solver is not available and licensed." ) @@ -4864,14 +4888,16 @@ def test_coefficient_matching_raises_error_4_3(self): # solve with PyROS dr_orders = [1, 2] for dr_order in dr_orders: - with self.assertRaisesRegex( + regex_assert_mgr = self.assertRaisesRegex( ValueError, expected_regex=( - "Equality constraint.*cannot be guaranteed to be robustly " - "feasible.*" + "Coefficient matching unsuccessful. See the solver logs." ), - ): - res = pyros_solver.solve( + ) + logging_intercept_mgr = LoggingIntercept(level=logging.ERROR) + + with regex_assert_mgr, logging_intercept_mgr as LOG: + pyros_solver.solve( model=m, first_stage_variables=[], second_stage_variables=[m.x1, m.x2, m.x3], @@ -4886,6 +4912,16 @@ def test_coefficient_matching_raises_error_4_3(self): robust_feasibility_tolerance=1e-4, ) + detailed_error_msg = LOG.getvalue() + self.assertRegex( + detailed_error_msg[:-1], + ( + r"Equality constraint.*cannot be guaranteed to " + r"be robustly feasible.*" + r"Consider editing this constraint.*" + ), + ) + def test_coefficient_matching_robust_infeasible_proof_in_pyros(self): # Write the deterministic Pyomo model m = ConcreteModel() @@ -5011,28 +5047,42 @@ def test_bypass_global_separation(self): global_subsolver = SolverFactory("baron") # Call the PyROS solver - results = pyros_solver.solve( - model=m, - first_stage_variables=[m.x1], - second_stage_variables=[m.x2], - uncertain_params=[m.u], - uncertainty_set=interval, - local_solver=local_subsolver, - global_solver=global_subsolver, - options={ - "objective_focus": ObjectiveType.worst_case, - "solve_master_globally": True, - "decision_rule_order": 0, - "bypass_global_separation": True, - }, - ) + with LoggingIntercept(level=logging.WARNING) as LOG: + results = pyros_solver.solve( + model=m, + first_stage_variables=[m.x1], + second_stage_variables=[m.x2], + uncertain_params=[m.u], + uncertainty_set=interval, + local_solver=local_subsolver, + global_solver=global_subsolver, + options={ + "objective_focus": ObjectiveType.worst_case, + "solve_master_globally": True, + "decision_rule_order": 0, + "bypass_global_separation": True, + }, + ) + # check termination robust optimal self.assertEqual( results.pyros_termination_condition, pyrosTerminationCondition.robust_optimal, msg="Returned termination condition is not return robust_optimal.", ) + # since robust optimal, we also expect warning-level logger + # message about bypassing of global separation subproblems + warning_msgs = LOG.getvalue() + self.assertRegex( + warning_msgs, + ( + r".*Option to bypass global separation was chosen\. " + r"Robust feasibility and optimality of the reported " + r"solution are not guaranteed\." + ), + ) + @unittest.skipUnless( baron_available and baron_license_is_valid, @@ -5199,12 +5249,26 @@ def test_multiple_objs(self): # and solve again m.obj_max = Objective(expr=-m.obj.expr, sense=pyo_max) m.obj.deactivate() - res = pyros_solver.solve(**solve_kwargs) + max_obj_res = pyros_solver.solve(**solve_kwargs) # check active objectives self.assertEqual(len(list(m.component_data_objects(Objective, active=True))), 1) self.assertTrue(m.obj_max.active) + self.assertTrue( + math.isclose( + res.final_objective_value, + -max_obj_res.final_objective_value, + abs_tol=2e-4, # 2x the default robust feasibility tolerance + ), + msg=( + f"Robust optimal objective value {res.final_objective_value} " + "for problem with minimization objective not close to " + f"negative of value {max_obj_res.final_objective_value} " + "of equivalent maximization objective." + ), + ) + class testModelIdentifyObjectives(unittest.TestCase): """ @@ -5631,5 +5695,420 @@ def test_two_stg_mod_with_intersection_set(self): ) +class TestIterationLogRecord(unittest.TestCase): + """ + Test the PyROS `IterationLogRecord` class. + """ + + def test_log_header(self): + """Test method for logging iteration log table header.""" + ans = ( + "------------------------------------------------------------------------------\n" + "Itn Objective 1-Stg Shift DR Shift #CViol Max Viol Wall Time (s)\n" + "------------------------------------------------------------------------------\n" + ) + with LoggingIntercept(level=logging.INFO) as LOG: + IterationLogRecord.log_header(logger.info) + + self.assertEqual( + LOG.getvalue(), + ans, + msg="Messages logged for iteration table header do not match expected result", + ) + + def test_log_standard_iter_record(self): + """Test logging function for PyROS IterationLogRecord.""" + + # for some fields, we choose floats with more than four + # four decimal points to ensure rounding also matches + iter_record = IterationLogRecord( + iteration=4, + objective=1.234567, + first_stage_var_shift=2.3456789e-8, + dr_var_shift=3.456789e-7, + num_violated_cons=10, + max_violation=7.654321e-3, + elapsed_time=21.2, + dr_polishing_success=True, + all_sep_problems_solved=True, + global_separation=False, + ) + + # now check record logged as expected + ans = ( + "4 1.2346e+00 2.3457e-08 3.4568e-07 10 7.6543e-03 " + "21.200 \n" + ) + with LoggingIntercept(level=logging.INFO) as LOG: + iter_record.log(logger.info) + result = LOG.getvalue() + + self.assertEqual( + ans, + result, + msg="Iteration log record message does not match expected result", + ) + + def test_log_iter_record_polishing_failed(self): + """Test iteration log record in event of polishing failure.""" + # for some fields, we choose floats with more than four + # four decimal points to ensure rounding also matches + iter_record = IterationLogRecord( + iteration=4, + objective=1.234567, + first_stage_var_shift=2.3456789e-8, + dr_var_shift=3.456789e-7, + num_violated_cons=10, + max_violation=7.654321e-3, + elapsed_time=21.2, + dr_polishing_success=False, + all_sep_problems_solved=True, + global_separation=False, + ) + + # now check record logged as expected + ans = ( + "4 1.2346e+00 2.3457e-08 3.4568e-07* 10 7.6543e-03 " + "21.200 \n" + ) + with LoggingIntercept(level=logging.INFO) as LOG: + iter_record.log(logger.info) + result = LOG.getvalue() + + self.assertEqual( + ans, + result, + msg="Iteration log record message does not match expected result", + ) + + def test_log_iter_record_global_separation(self): + """ + Test iteration log record in event global separation performed. + In this case, a 'g' should be appended to the max violation + reported. Useful in the event neither local nor global separation + was bypassed. + """ + # for some fields, we choose floats with more than four + # four decimal points to ensure rounding also matches + iter_record = IterationLogRecord( + iteration=4, + objective=1.234567, + first_stage_var_shift=2.3456789e-8, + dr_var_shift=3.456789e-7, + num_violated_cons=10, + max_violation=7.654321e-3, + elapsed_time=21.2, + dr_polishing_success=True, + all_sep_problems_solved=True, + global_separation=True, + ) + + # now check record logged as expected + ans = ( + "4 1.2346e+00 2.3457e-08 3.4568e-07 10 7.6543e-03g " + "21.200 \n" + ) + with LoggingIntercept(level=logging.INFO) as LOG: + iter_record.log(logger.info) + result = LOG.getvalue() + + self.assertEqual( + ans, + result, + msg="Iteration log record message does not match expected result", + ) + + def test_log_iter_record_not_all_sep_solved(self): + """ + Test iteration log record in event not all separation problems + were solved successfully. This may have occurred if the PyROS + solver time limit was reached, or the user-provides subordinate + optimizer(s) were unable to solve a separation subproblem + to an acceptable level. + A '+' should be appended to the number of performance constraints + found to be violated. + """ + # for some fields, we choose floats with more than four + # four decimal points to ensure rounding also matches + iter_record = IterationLogRecord( + iteration=4, + objective=1.234567, + first_stage_var_shift=2.3456789e-8, + dr_var_shift=3.456789e-7, + num_violated_cons=10, + max_violation=7.654321e-3, + elapsed_time=21.2, + dr_polishing_success=True, + all_sep_problems_solved=False, + global_separation=False, + ) + + # now check record logged as expected + ans = ( + "4 1.2346e+00 2.3457e-08 3.4568e-07 10+ 7.6543e-03 " + "21.200 \n" + ) + with LoggingIntercept(level=logging.INFO) as LOG: + iter_record.log(logger.info) + result = LOG.getvalue() + + self.assertEqual( + ans, + result, + msg="Iteration log record message does not match expected result", + ) + + def test_log_iter_record_all_special(self): + """ + Test iteration log record in event DR polishing and global + separation failed. + """ + # for some fields, we choose floats with more than four + # four decimal points to ensure rounding also matches + iter_record = IterationLogRecord( + iteration=4, + objective=1.234567, + first_stage_var_shift=2.3456789e-8, + dr_var_shift=3.456789e-7, + num_violated_cons=10, + max_violation=7.654321e-3, + elapsed_time=21.2, + dr_polishing_success=False, + all_sep_problems_solved=False, + global_separation=True, + ) + + # now check record logged as expected + ans = ( + "4 1.2346e+00 2.3457e-08 3.4568e-07* 10+ 7.6543e-03g " + "21.200 \n" + ) + with LoggingIntercept(level=logging.INFO) as LOG: + iter_record.log(logger.info) + result = LOG.getvalue() + + self.assertEqual( + ans, + result, + msg="Iteration log record message does not match expected result", + ) + + def test_log_iter_record_attrs_none(self): + """ + Test logging of iteration record in event some + attributes are of value `None`. In this case, a '-' + should be printed in lieu of a numerical value. + Example where this occurs: the first iteration, + in which there is no first-stage shift or DR shift. + """ + # for some fields, we choose floats with more than four + # four decimal points to ensure rounding also matches + iter_record = IterationLogRecord( + iteration=0, + objective=-1.234567, + first_stage_var_shift=None, + dr_var_shift=None, + num_violated_cons=10, + max_violation=7.654321e-3, + elapsed_time=21.2, + dr_polishing_success=True, + all_sep_problems_solved=False, + global_separation=True, + ) + + # now check record logged as expected + ans = ( + "0 -1.2346e+00 - - 10+ 7.6543e-03g " + "21.200 \n" + ) + with LoggingIntercept(level=logging.INFO) as LOG: + iter_record.log(logger.info) + result = LOG.getvalue() + + self.assertEqual( + ans, + result, + msg="Iteration log record message does not match expected result", + ) + + +class TestROSolveResults(unittest.TestCase): + """ + Test PyROS solver results object. + """ + + def test_ro_solve_results_str(self): + """ + Test string representation of RO solve results object. + """ + res = ROSolveResults( + config=SolverFactory("pyros").CONFIG(), + iterations=4, + final_objective_value=123.456789, + time=300.34567, + pyros_termination_condition=pyrosTerminationCondition.robust_optimal, + ) + ans = ( + "Termination stats:\n" + " Iterations : 4\n" + " Solve time (wall s) : 300.346\n" + " Final objective value : 1.2346e+02\n" + " Termination condition : pyrosTerminationCondition.robust_optimal" + ) + self.assertEqual( + str(res), + ans, + msg=( + "String representation of PyROS results object does not " + "match expected value" + ), + ) + + def test_ro_solve_results_str_attrs_none(self): + """ + Test string representation of PyROS solve results in event + one of the printed attributes is of value `None`. + This may occur at instantiation or, for example, + whenever the PyROS solver confirms robust infeasibility through + coefficient matching. + """ + res = ROSolveResults( + config=SolverFactory("pyros").CONFIG(), + iterations=0, + final_objective_value=None, + time=300.34567, + pyros_termination_condition=pyrosTerminationCondition.robust_optimal, + ) + ans = ( + "Termination stats:\n" + " Iterations : 0\n" + " Solve time (wall s) : 300.346\n" + " Final objective value : None\n" + " Termination condition : pyrosTerminationCondition.robust_optimal" + ) + self.assertEqual( + str(res), + ans, + msg=( + "String representation of PyROS results object does not " + "match expected value" + ), + ) + + +class TestPyROSSolverLogIntros(unittest.TestCase): + """ + Test logging of introductory information by PyROS solver. + """ + + def test_log_config(self): + """ + Test method for logging PyROS solver config dict. + """ + pyros_solver = SolverFactory("pyros") + config = pyros_solver.CONFIG(dict(nominal_uncertain_param_vals=[0.5])) + with LoggingIntercept(level=logging.INFO) as LOG: + pyros_solver._log_config(logger=logger, config=config, level=logging.INFO) + + ans = ( + "Solver options:\n" + " time_limit=None\n" + " keepfiles=False\n" + " tee=False\n" + " load_solution=True\n" + " objective_focus=\n" + " nominal_uncertain_param_vals=[0.5]\n" + " decision_rule_order=0\n" + " solve_master_globally=False\n" + " max_iter=-1\n" + " robust_feasibility_tolerance=0.0001\n" + " separation_priority_order={}\n" + " progress_logger=\n" + " backup_local_solvers=[]\n" + " backup_global_solvers=[]\n" + " subproblem_file_directory=None\n" + " bypass_local_separation=False\n" + " bypass_global_separation=False\n" + " p_robustness={}\n" + "-" * 78 + "\n" + ) + + logged_str = LOG.getvalue() + self.assertEqual( + logged_str, + ans, + msg=( + "Logger output for PyROS solver config (default case) " + "does not match expected result." + ), + ) + + def test_log_intro(self): + """ + Test logging of PyROS solver introductory messages. + """ + pyros_solver = SolverFactory("pyros") + with LoggingIntercept(level=logging.INFO) as LOG: + pyros_solver._log_intro(logger=logger, level=logging.INFO) + + intro_msgs = LOG.getvalue() + + # last character should be newline; disregard it + intro_msg_lines = intro_msgs.split("\n")[:-1] + + # check number of lines is as expected + self.assertEqual( + len(intro_msg_lines), + 14, + msg=( + "PyROS solver introductory message does not contain" + "the expected number of lines." + ), + ) + + # first and last lines of the introductory section + self.assertEqual(intro_msg_lines[0], "=" * 78) + self.assertEqual(intro_msg_lines[-1], "=" * 78) + + # check regex main text + self.assertRegex( + " ".join(intro_msg_lines[1:-1]), + r"PyROS: The Pyomo Robust Optimization Solver, v.* \(IDAES\)\.", + ) + + def test_log_disclaimer(self): + """ + Test logging of PyROS solver disclaimer messages. + """ + pyros_solver = SolverFactory("pyros") + with LoggingIntercept(level=logging.INFO) as LOG: + pyros_solver._log_disclaimer(logger=logger, level=logging.INFO) + + disclaimer_msgs = LOG.getvalue() + + # last character should be newline; disregard it + disclaimer_msg_lines = disclaimer_msgs.split("\n")[:-1] + + # check number of lines is as expected + self.assertEqual( + len(disclaimer_msg_lines), + 5, + msg=( + "PyROS solver disclaimer message does not contain" + "the expected number of lines." + ), + ) + + # regex first line of disclaimer section + self.assertRegex(disclaimer_msg_lines[0], r"=.* DISCLAIMER .*=") + # check last line of disclaimer section + self.assertEqual(disclaimer_msg_lines[-1], "=" * 78) + + # check regex main text + self.assertRegex( + " ".join(disclaimer_msg_lines[1:-1]), + r"PyROS is still under development.*ticket at.*", + ) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/pyros/util.py b/pyomo/contrib/pyros/util.py index 2c1a309ced3..a956e17b089 100644 --- a/pyomo/contrib/pyros/util.py +++ b/pyomo/contrib/pyros/util.py @@ -41,6 +41,8 @@ import logging from pprint import pprint import math +from pyomo.common.timing import HierarchicalTimer +from pyomo.common.log import Preformatted # Tolerances used in the code @@ -50,6 +52,135 @@ COEFF_MATCH_ABS_TOL = 0 ABS_CON_CHECK_FEAS_TOL = 1e-5 TIC_TOC_SOLVE_TIME_ATTR = "pyros_tic_toc_time" +DEFAULT_LOGGER_NAME = "pyomo.contrib.pyros" + + +class TimingData: + """ + PyROS solver timing data object. + + Implemented as a wrapper around `common.timing.HierarchicalTimer`, + with added functionality for enforcing a standardized + hierarchy of identifiers. + + Attributes + ---------- + hierarchical_timer_full_ids : set of str + (Class attribute.) Valid identifiers for use with + the encapsulated hierarchical timer. + """ + + hierarchical_timer_full_ids = { + "main", + "main.preprocessing", + "main.master_feasibility", + "main.master", + "main.dr_polishing", + "main.local_separation", + "main.global_separation", + } + + def __init__(self): + """Initialize self (see class docstring).""" + self._hierarchical_timer = HierarchicalTimer() + + def __str__(self): + """ + String representation of `self`. Currently + returns the string representation of `self.hierarchical_timer`. + + Returns + ------- + str + String representation. + """ + return self._hierarchical_timer.__str__() + + def _validate_full_identifier(self, full_identifier): + """ + Validate identifier for hierarchical timer. + + Parameters + ---------- + full_identifier : str + Identifier to validate. + + Raises + ------ + ValueError + If identifier not in `TimingData.hierarchical_timer_full_ids`. + """ + if full_identifier not in self.hierarchical_timer_full_ids: + raise ValueError( + "PyROS timing data object does not support timing ID: " + f"{full_identifier}." + ) + + def start_timer(self, full_identifier): + """ + Start timer for `self.hierarchical_timer`. + + Parameters + ---------- + full_identifier : str + Full identifier for the timer to be started. + Must be an entry of + `TimingData.hierarchical_timer_full_ids`. + """ + self._validate_full_identifier(full_identifier) + identifier = full_identifier.split(".")[-1] + return self._hierarchical_timer.start(identifier=identifier) + + def stop_timer(self, full_identifier): + """ + Stop timer for `self.hierarchical_timer`. + + Parameters + ---------- + full_identifier : str + Full identifier for the timer to be stopped. + Must be an entry of + `TimingData.hierarchical_timer_full_ids`. + """ + self._validate_full_identifier(full_identifier) + identifier = full_identifier.split(".")[-1] + return self._hierarchical_timer.stop(identifier=identifier) + + def get_total_time(self, full_identifier): + """ + Get total time spent with identifier active. + + Parameters + ---------- + full_identifier : str + Full identifier for the timer of interest. + + Returns + ------- + float + Total time spent with identifier active. + """ + return self._hierarchical_timer.get_total_time(identifier=full_identifier) + + def get_main_elapsed_time(self): + """ + Get total time elapsed for main timer of + the HierarchicalTimer contained in self. + + Returns + ------- + float + Total elapsed time. + + Note + ---- + This method is meant for use while the main timer is active. + Otherwise, use ``self.get_total_time("main")``. + """ + # clean? + return self._hierarchical_timer.timers["main"].tic_toc.toc( + msg=None, delta=False + ) '''Code borrowed from gdpopt: time_code, get_main_elapsed_time, a_logger.''' @@ -57,31 +188,33 @@ @contextmanager def time_code(timing_data_obj, code_block_name, is_main_timer=False): - """Starts timer at entry, stores elapsed time at exit + """ + Starts timer at entry, stores elapsed time at exit. + + Parameters + ---------- + timing_data_obj : TimingData + Timing data object. + code_block_name : str + Name of code block being timed. If `is_main_timer=True`, the start time is stored in the timing_data_obj, allowing calculation of total elapsed time 'on the fly' (e.g. to enforce a time limit) using `get_main_elapsed_time(timing_data_obj)`. """ + # initialize tic toc timer + timing_data_obj.start_timer(code_block_name) + start_time = timeit.default_timer() if is_main_timer: timing_data_obj.main_timer_start_time = start_time yield - elapsed_time = timeit.default_timer() - start_time - prev_time = timing_data_obj.get(code_block_name, 0) - timing_data_obj[code_block_name] = prev_time + elapsed_time + timing_data_obj.stop_timer(code_block_name) def get_main_elapsed_time(timing_data_obj): """Returns the time since entering the main `time_code` context""" - current_time = timeit.default_timer() - try: - return current_time - timing_data_obj.main_timer_start_time - except AttributeError as e: - if 'main_timer_start_time' in str(e): - raise AttributeError( - "You need to be in a 'time_code' context to use `get_main_elapsed_time()`." - ) + return timing_data_obj.get_main_elapsed_time() def adjust_solver_time_settings(timing_data_obj, solver, config): @@ -223,10 +356,107 @@ def revert_solver_max_time_adjustment( del solver.options[options_key] +class PreformattedLogger(logging.Logger): + """ + A specialized logger object designed to cast log messages + to Pyomo `Preformatted` objects prior to logging the messages. + Useful for circumventing the formatters of the standard Pyomo + logger in the event an instance is a descendant of the Pyomo + logger. + """ + + def critical(self, msg, *args, **kwargs): + """ + Preformat and log ``msg % args`` with severity + `logging.CRITICAL`. + """ + return super(PreformattedLogger, self).critical( + Preformatted(msg % args if args else msg), **kwargs + ) + + def error(self, msg, *args, **kwargs): + """ + Preformat and log ``msg % args`` with severity + `logging.ERROR`. + """ + return super(PreformattedLogger, self).error( + Preformatted(msg % args if args else msg), **kwargs + ) + + def warning(self, msg, *args, **kwargs): + """ + Preformat and log ``msg % args`` with severity + `logging.WARNING`. + """ + return super(PreformattedLogger, self).warning( + Preformatted(msg % args if args else msg), **kwargs + ) + + def info(self, msg, *args, **kwargs): + """ + Preformat and log ``msg % args`` with severity + `logging.INFO`. + """ + return super(PreformattedLogger, self).info( + Preformatted(msg % args if args else msg), **kwargs + ) + + def debug(self, msg, *args, **kwargs): + """ + Preformat and log ``msg % args`` with severity + `logging.DEBUG`. + """ + return super(PreformattedLogger, self).debug( + Preformatted(msg % args if args else msg), **kwargs + ) + + def log(self, level, msg, *args, **kwargs): + """ + Preformat and log ``msg % args`` with integer + severity `level`. + """ + return super(PreformattedLogger, self).log( + level, Preformatted(msg % args if args else msg), **kwargs + ) + + +def setup_pyros_logger(name=DEFAULT_LOGGER_NAME): + """ + Set up pyros logger. + """ + # default logger: INFO level, with preformatted messages + current_logger_class = logging.getLoggerClass() + logging.setLoggerClass(PreformattedLogger) + logger = logging.getLogger(DEFAULT_LOGGER_NAME) + logger.setLevel(logging.INFO) + logging.setLoggerClass(current_logger_class) + + return logger + + def a_logger(str_or_logger): - """Returns a logger when passed either a logger name or logger object.""" + """ + Standardize a string or logger object to a logger object. + + Parameters + ---------- + str_or_logger : str or logging.Logger + String or logger object to normalize. + + Returns + ------- + logging.Logger + If `str_or_logger` is of type `logging.Logger`,then + `str_or_logger` is returned. + Otherwise, ``logging.getLogger(str_or_logger)`` + is returned. In the event `str_or_logger` is + the name of the default PyROS logger, the logger level + is set to `logging.INFO`, and a `PreformattedLogger` + instance is returned in lieu of a standard `Logger` + instance. + """ if isinstance(str_or_logger, logging.Logger): - return str_or_logger + return logging.getLogger(str_or_logger.name) else: return logging.getLogger(str_or_logger) @@ -271,6 +501,25 @@ class pyrosTerminationCondition(Enum): time_out = 5 """Maximum allowable time exceeded.""" + @property + def message(self): + """ + str : Message associated with a given PyROS + termination condition. + """ + message_dict = { + self.robust_optimal: "Robust optimal solution identified.", + self.robust_feasible: "Robust feasible solution identified.", + self.robust_infeasible: "Problem is robust infeasible.", + self.time_out: "Maximum allowable time exceeded.", + self.max_iter: "Maximum number of iterations reached.", + self.subsolver_error: ( + "Subordinate optimizer(s) could not solve a subproblem " + "to an acceptable status." + ), + } + return message_dict[self] + class SeparationStrategy(Enum): all_violations = auto() @@ -1383,111 +1632,193 @@ def process_termination_condition_master_problem(config, results): ) -def output_logger(config, **kwargs): - ''' - All user returned messages (termination conditions, runtime errors) are here - Includes when - "sub-solver %s returned status infeasible..." - :return: - ''' +class IterationLogRecord: + """ + PyROS solver iteration log record. - # === PREAMBLE + LICENSING - # Version printing - if "preamble" in kwargs: - if kwargs["preamble"]: - version = str(kwargs["version"]) - preamble = ( - "===========================================================================================\n" - "PyROS: Pyomo Robust Optimization Solver v.%s \n" - "Developed by: Natalie M. Isenberg (1), Jason A. F. Sherman (1), \n" - " John D. Siirola (2), Chrysanthos E. Gounaris (1) \n" - "(1) Carnegie Mellon University, Department of Chemical Engineering \n" - "(2) Sandia National Laboratories, Center for Computing Research\n\n" - "The developers gratefully acknowledge support from the U.S. Department of Energy's \n" - "Institute for the Design of Advanced Energy Systems (IDAES) \n" - "===========================================================================================" - % version - ) - print(preamble) - # === DISCLAIMER - if "disclaimer" in kwargs: - if kwargs["disclaimer"]: - print( - "======================================== DISCLAIMER =======================================\n" - "PyROS is still under development. \n" - "Please provide feedback and/or report any issues by opening a Pyomo ticket.\n" - "===========================================================================================\n" - ) - # === ALL LOGGER RETURN MESSAGES - if "bypass_global_separation" in kwargs: - if kwargs["bypass_global_separation"]: - config.progress_logger.info( - "NOTE: Option to bypass global separation was chosen. " - "Robust feasibility and optimality of the reported " - "solution are not guaranteed." - ) - if "robust_optimal" in kwargs: - if kwargs["robust_optimal"]: - config.progress_logger.info( - 'Robust optimal solution identified. Exiting PyROS.' - ) + Parameters + ---------- + iteration : int or None, optional + Iteration number. + objective : int or None, optional + Master problem objective value. + Note: if the sense of the original model is maximization, + then this is the negative of the objective value + of the original model. + first_stage_var_shift : float or None, optional + Infinity norm of the difference between first-stage + variable vectors for the current and previous iterations. + dr_var_shift : float or None, optional + Infinity norm of the difference between decision rule + variable vectors for the current and previous iterations. + dr_polishing_success : bool or None, optional + True if DR polishing solved successfully, False otherwise. + num_violated_cons : int or None, optional + Number of performance constraints found to be violated + during separation step. + all_sep_problems_solved : int or None, optional + True if all separation problems were solved successfully, + False otherwise (such as if there was a time out, subsolver + error, or only a subset of the problems were solved due to + custom constraint prioritization). + global_separation : bool, optional + True if separation problems were solved with the subordinate + global optimizer(s), False otherwise. + max_violation : int or None + Maximum scaled violation of any performance constraint + found during separation step. + elapsed_time : float, optional + Total time elapsed up to the current iteration, in seconds. + + Attributes + ---------- + iteration : int or None + Iteration number. + objective : int or None + Master problem objective value. + Note: if the sense of the original model is maximization, + then this is the negative of the objective value + of the original model. + first_stage_var_shift : float or None + Infinity norm of the difference between first-stage + variable vectors for the current and previous iterations. + dr_var_shift : float or None + Infinity norm of the difference between decision rule + variable vectors for the current and previous iterations. + dr_polishing_success : bool or None + True if DR polishing solved successfully, False otherwise. + num_violated_cons : int or None + Number of performance constraints found to be violated + during separation step. + all_sep_problems_solved : int or None + True if all separation problems were solved successfully, + False otherwise (such as if there was a time out, subsolver + error, or only a subset of the problems were solved due to + custom constraint prioritization). + global_separation : bool + True if separation problems were solved with the subordinate + global optimizer(s), False otherwise. + max_violation : int or None + Maximum scaled violation of any performance constraint + found during separation step. + elapsed_time : float + Total time elapsed up to the current iteration, in seconds. + """ - if "robust_feasible" in kwargs: - if kwargs["robust_feasible"]: - config.progress_logger.info( - 'Robust feasible solution identified. Exiting PyROS.' - ) + _LINE_LENGTH = 78 + _ATTR_FORMAT_LENGTHS = { + "iteration": 5, + "objective": 13, + "first_stage_var_shift": 13, + "dr_var_shift": 13, + "num_violated_cons": 8, + "max_violation": 13, + "elapsed_time": 13, + } + _ATTR_HEADER_NAMES = { + "iteration": "Itn", + "objective": "Objective", + "first_stage_var_shift": "1-Stg Shift", + "dr_var_shift": "DR Shift", + "num_violated_cons": "#CViol", + "max_violation": "Max Viol", + "elapsed_time": "Wall Time (s)", + } + + def __init__( + self, + iteration, + objective, + first_stage_var_shift, + dr_var_shift, + dr_polishing_success, + num_violated_cons, + all_sep_problems_solved, + global_separation, + max_violation, + elapsed_time, + ): + """Initialize self (see class docstring).""" + self.iteration = iteration + self.objective = objective + self.first_stage_var_shift = first_stage_var_shift + self.dr_var_shift = dr_var_shift + self.dr_polishing_success = dr_polishing_success + self.num_violated_cons = num_violated_cons + self.all_sep_problems_solved = all_sep_problems_solved + self.global_separation = global_separation + self.max_violation = max_violation + self.elapsed_time = elapsed_time + + def get_log_str(self): + """Get iteration log string.""" + attrs = [ + "iteration", + "objective", + "first_stage_var_shift", + "dr_var_shift", + "num_violated_cons", + "max_violation", + "elapsed_time", + ] + return "".join(self._format_record_attr(attr) for attr in attrs) + + def _format_record_attr(self, attr_name): + """Format attribute record for logging.""" + attr_val = getattr(self, attr_name) + if attr_val is None: + fmt_str = f"<{self._ATTR_FORMAT_LENGTHS[attr_name]}s" + return f"{'-':{fmt_str}}" + else: + attr_val_fstrs = { + "iteration": "f'{attr_val:d}'", + "objective": "f'{attr_val: .4e}'", + "first_stage_var_shift": "f'{attr_val:.4e}'", + "dr_var_shift": "f'{attr_val:.4e}'", + "num_violated_cons": "f'{attr_val:d}'", + "max_violation": "f'{attr_val:.4e}'", + "elapsed_time": "f'{attr_val:.3f}'", + } + + # qualifier for DR polishing and separation columns + if attr_name == "dr_var_shift": + qual = "*" if not self.dr_polishing_success else "" + elif attr_name == "num_violated_cons": + qual = "+" if not self.all_sep_problems_solved else "" + elif attr_name == "max_violation": + qual = "g" if self.global_separation else "" + else: + qual = "" - if "robust_infeasible" in kwargs: - if kwargs["robust_infeasible"]: - config.progress_logger.info('Robust infeasible problem. Exiting PyROS.') - - if "time_out" in kwargs: - if kwargs["time_out"]: - config.progress_logger.info( - 'PyROS was unable to identify robust solution ' - 'before exceeding time limit of %s seconds. ' - 'Consider increasing the time limit via option time_limit.' - % config.time_limit - ) + attr_val_str = f"{eval(attr_val_fstrs[attr_name])}{qual}" - if "max_iter" in kwargs: - if kwargs["max_iter"]: - config.progress_logger.info( - 'PyROS was unable to identify robust solution ' - 'within %s iterations of the GRCS algorithm. ' - 'Consider increasing the iteration limit via option max_iter.' - % config.max_iter - ) + return f"{attr_val_str:{f'<{self._ATTR_FORMAT_LENGTHS[attr_name]}'}}" - if "master_error" in kwargs: - if kwargs["master_error"]: - status_dict = kwargs["status_dict"] - filename = kwargs["filename"] # solver name to solver termination condition - if kwargs["iteration"] == 0: - raise AttributeError( - "User-supplied solver(s) could not solve the deterministic model. " - "Returned termination conditions were: %s" - "Please ensure deterministic model is solvable by at least one of the supplied solvers. " - "Exiting PyROS." % pprint(status_dict, width=1) - ) - config.progress_logger.info( - "User-supplied solver(s) could not solve the master model at iteration %s.\n" - "Returned termination conditions were: %s\n" - "For debugging, this problem has been written to a GAMS file titled %s. Exiting PyROS." - % (kwargs["iteration"], pprint(status_dict), filename) - ) - if "separation_error" in kwargs: - if kwargs["separation_error"]: - status_dict = kwargs["status_dict"] - filename = kwargs["filename"] - iteration = kwargs["iteration"] - obj = kwargs["objective"] - config.progress_logger.info( - "User-supplied solver(s) could not solve the separation problem at iteration %s under separation objective %s.\n" - "Returned termination conditions were: %s\n" - "For debugging, this problem has been written to a GAMS file titled %s. Exiting PyROS." - % (iteration, obj, pprint(status_dict, width=1), filename) - ) + def log(self, log_func, **log_func_kwargs): + """Log self.""" + log_str = self.get_log_str() + log_func(log_str, **log_func_kwargs) - return + @staticmethod + def get_log_header_str(): + """Get string for iteration log header.""" + fmt_lengths_dict = IterationLogRecord._ATTR_FORMAT_LENGTHS + header_names_dict = IterationLogRecord._ATTR_HEADER_NAMES + return "".join( + f"{header_names_dict[attr]:<{fmt_lengths_dict[attr]}s}" + for attr in fmt_lengths_dict + ) + + @staticmethod + def log_header(log_func, with_rules=True, **log_func_kwargs): + """Log header.""" + if with_rules: + IterationLogRecord.log_header_rule(log_func, **log_func_kwargs) + log_func(IterationLogRecord.get_log_header_str(), **log_func_kwargs) + if with_rules: + IterationLogRecord.log_header_rule(log_func, **log_func_kwargs) + + @staticmethod + def log_header_rule(log_func, fillchar="-", **log_func_kwargs): + """Log header rule.""" + log_func(fillchar * IterationLogRecord._LINE_LENGTH, **log_func_kwargs)