diff --git a/column_testing.py b/column_testing.py new file mode 100644 index 0000000..23dc6f8 --- /dev/null +++ b/column_testing.py @@ -0,0 +1,98 @@ +from __future__ import division + +# Import various modules and functions needed for the script +import csv # To handle CSV files +import logging # To keep logs for tracking +import os # To access the OS functionalities for file and directory handling +from math import ceil, fabs # Importing ceil and fabs functions from math +import time # To keep track of time taken to solve the optimization problem +import pyomo.environ as pe # To create and solve optimization models + +# Importing specific classes and functions from pyomo.environ +from pyomo.environ import ( + Block, + BooleanVar, + ConcreteModel, + Constraint, + NonNegativeReals, + Objective, + Param, + RangeSet, + Set, + SolverFactory, + Suffix, + TransformationFactory, + Var, + exactly, + land, + log, + lor, + minimize, + value, +) + +# Importing Disjunct and Disjunction classes from pyomo.gdp for creating generalized disjunctive programming models +from pyomo.gdp import Disjunct, Disjunction + +# Importing utility function to log infeasible constraints +from pyomo.util.infeasible import log_infeasible_constraints + +# Importing build_column function from gdp.column.gdp_column +from gdp.column.gdp_column import build_column + +if __name__ == "__main__": + # This part of the script initializes the variables that are going to be used throughout the script. + # NT: Number of trays in the distillation column + # timelimit: time limit for solving the optimization problem + # model_args: dictionary containing arguments for the distillation column model + # starting_point: List containing initial guesses for the number of trays and reflux. + # globaltee: Boolean indicating whether to display solver output. + # logging: Configuring the logging level to ERROR. This will avoid printing out warning messages. + + NT = 17 + time_limit = 900 # [s] + model_args = {'min_trays': 8, 'max_trays': NT, 'xD': 0.95, 'xB': 0.95} + starting_point = [NT - 2, 1] + globaltee = True + logging.basicConfig(level=logging.ERROR) + + nlps = ['knitro', 'baron'] + ks = ['Linf', 'L2'] + + m = build_column(**model_args) + + results = [] + + for solver in nlps: + for k in ks: + new_results={} + start_time = time.time() + result = pe.SolverFactory('gdpopt.ldsda').solve( + m, + minlp_solver='gams', + nlp_solver=solver, + tee=True, + starting_point=starting_point, + logical_constraint_list=[m.one_reflux, + m.one_boilup + ], + direction_norm=k, + time_limit=time_limit, + ) + end_time = time.time() + elapsed_time = end_time - start_time + new_results={"Solver": solver, "k": str(k), "Objective": pe.value(m.obj), "Time": elapsed_time} + results.append(new_results) + print(new_results) + + # Write the results to a CSV file + csv_file_path = 'column_testing.csv' + with open(csv_file_path, mode='w', newline='') as csv_file: + fieldnames = ['Solver', 'k', 'Objective', 'Time'] + writer = csv.DictWriter(csv_file, fieldnames=fieldnames) + + writer.writeheader() + for result in results: + writer.writerow(result) + + print(f"Results have been written to {csv_file_path}") \ No newline at end of file diff --git a/cstr_ldsda.csv b/cstr_ldsda.csv new file mode 100644 index 0000000..ef9315d --- /dev/null +++ b/cstr_ldsda.csv @@ -0,0 +1,105 @@ +NT,NLP_Solver,k,Objective,Time (s) +5,knitro,L2,3.130198132329715,3.7278621196746826 +5,baron,L2,3.130198132329715,2.775907039642334 +5,knitro,Linf,3.0620145765739197,3.7131240367889404 +5,baron,Linf,3.0620145765739197,3.63629412651062 +6,knitro,L2,3.130198132329715,5.418642520904541 +6,baron,L2,3.130198132329715,3.9788856506347656 +6,knitro,Linf,3.0072779739861715,4.427759170532227 +6,baron,Linf,3.0072779739861715,4.7388246059417725 +7,knitro,L2,3.130198132329715,4.075173616409302 +7,baron,L2,3.130198132329715,3.9738059043884277 +7,knitro,Linf,2.966047072054598,5.554377555847168 +7,baron,Linf,2.966047072054598,5.237062692642212 +8,knitro,L2,3.130198132329715,4.122750997543335 +8,baron,L2,3.130198132329715,3.7697317600250244 +8,knitro,Linf,2.9344051030872005,5.333848476409912 +8,baron,Linf,2.9344051030872005,5.706854581832886 +9,knitro,L2,3.130198132329715,3.3213605880737305 +9,baron,L2,3.130198132329715,4.077863931655884 +9,knitro,Linf,2.909533261913244,6.7725303173065186 +9,baron,Linf,2.909533261913244,7.81814432144165 +10,knitro,L2,3.130198132329715,4.754985809326172 +10,baron,L2,3.130198132329715,4.5682618618011475 +10,knitro,Linf,2.8895369685436294,8.39240550994873 +10,baron,Linf,2.8895369685436294,9.142786026000977 +11,knitro,L2,3.130198132329715,5.86026406288147 +11,baron,L2,3.130198132329715,5.832338094711304 +11,knitro,Linf,2.873140226203434,11.03695797920227 +11,baron,Linf,2.873140226203434,8.022186994552612 +12,knitro,L2,3.130198132329715,4.204470634460449 +12,baron,L2,3.130198132329715,4.034653186798096 +12,knitro,Linf,2.8594656076358045,10.365803003311157 +12,baron,Linf,2.8594656076358045,9.95569634437561 +13,knitro,L2,3.130198132329715,4.338276386260986 +13,baron,L2,3.130198132329715,3.8011109828948975 +13,knitro,Linf,2.8478946319864953,10.701895475387573 +13,baron,Linf,2.8478946319864953,9.914148092269897 +14,knitro,L2,3.130198132329715,4.2346556186676025 +14,baron,L2,3.130198132329715,5.409814834594727 +14,knitro,Linf,2.837980692760724,10.016036033630371 +14,baron,Linf,2.837980692760724,9.603132963180542 +15,knitro,L2,3.130198132329715,5.071369409561157 +15,baron,L2,3.130198132329715,5.214901924133301 +15,knitro,Linf,2.8293940630587144,12.931170463562012 +15,baron,Linf,2.8293940630587144,12.284100770950317 +16,knitro,L2,3.130198132329715,4.191612958908081 +16,baron,L2,3.130198132329715,5.170598983764648 +16,knitro,Linf,2.8218864094767833,15.052702188491821 +16,baron,Linf,2.8218864094767833,15.335492134094238 +17,knitro,L2,3.130198132329715,5.9339189529418945 +17,baron,L2,3.130198132329715,6.086316823959351 +17,knitro,Linf,2.8152673375976867,17.322603940963745 +17,baron,Linf,2.8152673375976867,17.62444019317627 +18,knitro,L2,3.130198132329715,6.724985122680664 +18,baron,L2,3.130198132329715,6.083144903182983 +18,knitro,Linf,2.8093885376330894,18.26237154006958 +18,baron,Linf,2.8093885376330894,16.826720476150513 +19,knitro,L2,3.130198132329715,8.952149868011475 +19,baron,L2,3.130198132329715,9.904943704605103 +19,knitro,Linf,2.80413282638558,27.2341411113739 +19,baron,Linf,2.80413282638558,23.826519012451172 +20,knitro,L2,3.130198132329715,8.070159673690796 +20,baron,L2,3.130198132329715,6.494313955307007 +20,knitro,Linf,2.79940642501676,17.605342864990234 +20,baron,Linf,2.79940642501676,16.106475353240967 +21,knitro,L2,3.130198132329715,4.7779927253723145 +21,baron,L2,3.130198132329715,5.472589492797852 +21,knitro,Linf,2.7951334164061983,17.953311681747437 +21,baron,Linf,2.7951334164061983,17.465514183044434 +22,knitro,L2,3.130198132329715,5.871582746505737 +22,baron,L2,3.130198132329715,7.571946144104004 +22,knitro,Linf,2.7912516990272183,19.04247546195984 +22,baron,Linf,2.7912516990272183,27.17935585975647 +23,knitro,L2,3.130198132329715,5.790415525436401 +23,baron,L2,3.130198132329715,5.989009141921997 +23,knitro,Linf,2.7877099905261713,22.043853759765625 +23,baron,Linf,2.7877099905261713,26.059494495391846 +24,knitro,L2,3.130198132329715,7.889803409576416 +24,baron,L2,3.130198132329715,7.548046588897705 +24,knitro,Linf,2.7844655755962724,45.06238317489624 +24,baron,Linf,2.7844655755962724,27.82550621032715 +25,knitro,L2,3.130198132329715,9.639538764953613 +25,baron,L2,3.130198132329715,7.886870861053467 +25,knitro,Linf,2.7814825926754,31.294180870056152 +25,baron,Linf,2.7814825926754,26.09557056427002 +26,knitro,L2,3.130198132329715,5.957033395767212 +26,baron,L2,3.130198132329715,5.945906400680542 +26,knitro,Linf,2.7787307151276415,31.193901777267456 +26,baron,Linf,2.7787307151276415,40.398841381073 +27,knitro,L2,3.130198132329715,14.585143089294434 +27,baron,L2,3.130198132329715,8.135599136352539 +27,knitro,Linf,2.7761841248417185,39.69969987869263 +27,baron,Linf,2.7761841248417185,38.46226906776428 +28,knitro,L2,3.130198132329715,8.921888828277588 +28,baron,L2,3.130198132329715,9.568353176116943 +28,knitro,Linf,2.773820706017508,43.28029990196228 +28,baron,Linf,2.773820706017508,55.194143772125244 +29,knitro,L2,3.130198132329715,8.783891439437866 +29,baron,L2,3.130198132329715,8.956895112991333 +29,knitro,Linf,2.7716214052136205,44.930012226104736 +29,baron,Linf,2.7716214052136205,40.95608878135681 +30,knitro,L2,3.130198132329715,9.60117769241333 +30,baron,L2,3.130198132329715,12.228340864181519 +30,knitro,Linf,2.7695697198462006,42.77383255958557 +30,baron,Linf,2.7695697198462006,38.65700149536133 diff --git a/cstr_testing.py b/cstr_testing.py new file mode 100644 index 0000000..2cb9972 --- /dev/null +++ b/cstr_testing.py @@ -0,0 +1,99 @@ +import csv +import logging +import os +from math import ceil, fabs +import time + +import pyomo.environ as pe +from pyomo.environ import SolverFactory, Suffix, value +from pyomo.gdp import Disjunct, Disjunction +from gdp.cstr.gdp_reactor import build_cstrs + +if __name__ == "__main__": + + # Results + NTs = range(5, 31, 1) + # NTs = [5] + time_limit = 900 + nlps = ['knitro', 'baron'] + starting_point = [1, 1] + + globaltee = True + + + transformations = ['bigm', 'hull'] + ks = ['L2', 'Linf'] + strategies = ['LOA', 'GLOA', 'LBB'] + + # Opening the CSV file here + with open('cstr_ldsda.csv', mode='w', newline='') as file: + writer = csv.writer(file) + # Write headers + writer.writerow(['NT', 'NLP_Solver', 'k', 'Objective', 'Time (s)']) + + # LD-SDA + for NT in NTs: + for k in ks: + for solver in nlps: + m = build_cstrs(NT) + new_results = {} + start_time = time.time() + result = pe.SolverFactory('gdpopt.ldsda').solve( + m, + minlp_solver='gams', + nlp_solver=solver, + tee=True, + starting_point=starting_point, + logical_constraint_list=[m.one_unreacted_feed, + m.one_recycle + ], + direction_norm=k, + time_limit=time_limit, + ) + + end_time = time.time() # End timing + elapsed_time = end_time - start_time # Calculate elapsed time + new_results = { + 'NT': NT, + 'Solver': solver, + 'k': str(k), + 'Objective': pe.value(m.obj), + 'Time': elapsed_time + } + writer.writerow([new_results['NT'], new_results['Solver'], new_results['k'], new_results['Objective'], new_results['Time']]) + print(new_results) + + # # D-SDA + # m = build_cstrs(NT) + # ext_ref = {m.YF: m.N, m.YR: m.N} + # get_external_information(m, ext_ref, tee=False) + + # for solver in nlps: + # for k in ks: + # for transformation in ['hull','bigm']: + # new_result = {} + # m_solved, _, _ = solve_with_dsda( + # model_function=build_cstrs, + # model_args={'NT': NT}, + # starting_point=starting_point, + # ext_dict=ext_ref, + # ext_logic=problem_logic_cstr, + # mip_transformation=True, + # transformation=transformation, + # k=k, + # provide_starting_initialization=True, + # feasible_model='cstr_' + str(NT), + # subproblem_solver=solver, + # subproblem_solver_options=nlp_opts[solver], + # iter_timelimit=timelimit, + # timelimit=timelimit, + # gams_output=False, + # tee=False, + # global_tee=False, + # ) + # new_result = {'Method': str('D-SDA_MIP_'+transformation), 'Approach': str('k='+k), 'Solver': solver, 'Objective': pe.value( + # m_solved.obj), 'Time': m_solved.dsda_time, 'Status': m_solved.dsda_status, 'User_time': m_solved.dsda_usertime, 'NT': NT} + # dict_data.append(new_result) + # print(new_result) + + diff --git a/gdp/column/column.py b/gdp/column/column.py index aabdeca..0873fad 100644 --- a/gdp/column/column.py +++ b/gdp/column/column.py @@ -28,7 +28,7 @@ from pyomo.common.errors import InfeasibleConstraintException from pyomo.contrib.fbbt.fbbt import fbbt -from pyomo.contrib.gdpopt.data_class import MasterProblemResult +#from pyomo.contrib.gdpopt.data_class import MasterProblemResult from pyomo.environ import ( Block, BooleanVar, diff --git a/gdp/column/gdp_column.py b/gdp/column/gdp_column.py index 2d253db..65f282f 100644 --- a/gdp/column/gdp_column.py +++ b/gdp/column/gdp_column.py @@ -21,7 +21,7 @@ # Imports from the Pyomo library for building and solving optimization problems. from pyomo.common.errors import InfeasibleConstraintException from pyomo.contrib.fbbt.fbbt import fbbt -from pyomo.contrib.gdpopt.data_class import MasterProblemResult +#from pyomo.contrib.gdpopt.data_class import MasterProblemResult from pyomo.core.base.misc import display from pyomo.core.plugins.transform.logical_to_linear import ( update_boolean_vars_from_binary, diff --git a/gdp/dsda/dsda_functions.py b/gdp/dsda/dsda_functions.py index b4e3015..05182cd 100644 --- a/gdp/dsda/dsda_functions.py +++ b/gdp/dsda/dsda_functions.py @@ -11,7 +11,6 @@ from gdp.dsda.model_serializer import StoreSpec, from_json, to_json from pyomo.common.errors import InfeasibleConstraintException from pyomo.contrib.fbbt.fbbt import fbbt -from pyomo.contrib.gdpopt.data_class import MasterProblemResult from pyomo.core.base.misc import display from pyomo.core.plugins.transform.logical_to_linear import \ update_boolean_vars_from_binary diff --git a/main_column.py b/main_column.py index 5aaa317..25bcda0 100644 --- a/main_column.py +++ b/main_column.py @@ -392,7 +392,7 @@ def problem_logic_column(m): model_args=model_args, starting_point=starting_point, ext_dict=ext_ref, - mip_transformation=True, + mip_transformation=False, transformation=transformation, ext_logic=problem_logic_column, k=k, diff --git a/main_cstr.py b/main_cstr.py index d5682f5..af0aad5 100644 --- a/main_cstr.py +++ b/main_cstr.py @@ -255,7 +255,7 @@ def problem_logic_cstr(m): starting_point=starting_point, ext_dict=ext_ref, ext_logic=problem_logic_cstr, - mip_transformation=True, + mip_transformation=False, transformation=transformation, k=k, provide_starting_initialization=True, diff --git a/main_small_batch.py b/main_small_batch.py index 19604eb..8cb1a83 100644 --- a/main_small_batch.py +++ b/main_small_batch.py @@ -152,7 +152,7 @@ def problem_logic_batch(m): model_args={}, starting_point=starting_point, ext_dict=ext_ref, - mip_transformation=True, + mip_transformation=False, transformation=transformation, ext_logic=problem_logic_batch, k=k, @@ -163,7 +163,7 @@ def problem_logic_batch(m): iter_timelimit=timelimit, timelimit=timelimit, gams_output=False, - tee=globaltee, + tee=False, global_tee=globaltee, ) new_result = {'Method': str('D-SDA_MIP_'+transformation), 'Approach': str('k='+k), 'Solver': solver, 'Objective': pe.value( diff --git a/results.csv b/results.csv new file mode 100644 index 0000000..fb865f4 --- /dev/null +++ b/results.csv @@ -0,0 +1,5 @@ +Solver,k,Objective,Time +knitro,Linf,19346.124801171078,76.94611406326294 +knitro,L2,19449.85151196701,13.004862308502197 +baron,Linf,19346.124802225793,18.10545325279236 +baron,L2,19449.851511967005,13.560934066772461 diff --git a/small_batch_ldsda.csv b/small_batch_ldsda.csv new file mode 100644 index 0000000..3dd449d --- /dev/null +++ b/small_batch_ldsda.csv @@ -0,0 +1,5 @@ +NLP_Solver,k,Objective,Time (s) +knitro,L2,167427.65710272908,4.435163736343384 +knitro,Linf,167427.65710272905,11.571666479110718 +baron,L2,167427.65710272905,5.39567756652832 +baron,Linf,167427.65710272902,11.289292573928833 diff --git a/small_batch_testing.py b/small_batch_testing.py new file mode 100644 index 0000000..b968ac3 --- /dev/null +++ b/small_batch_testing.py @@ -0,0 +1,64 @@ +import csv +import logging +import os +from math import ceil, fabs +import time + +import pyomo.environ as pe +from pyomo.gdp import Disjunct, Disjunction +from pyomo.util.infeasible import log_infeasible_constraints + +from gdp.small_batch.gdp_small_batch import build_small_batch + +if __name__ == "__main__": + # Inputs + time_limit = 900 + model_args = {} + starting_point = [3, 3, 3] + + globaltee = True + + nlps = ['knitro', 'baron'] + + transformations = ['bigm', 'hull'] + ks = ['L2', 'Linf'] + + m = build_small_batch() + + with open('small_batch_ldsda.csv', mode='w', newline='') as file: + writer = csv.writer(file) + # Write headers + writer.writerow(['NLP_Solver', 'k', 'Objective', 'Time (s)']) + + for solver in nlps: + for k in ks: + new_results={} + start_time = time.time() + result = pe.SolverFactory('gdpopt.ldsda').solve( + m, + minlp_solver='gams', + nlp_solver=solver, + tee=True, + starting_point=[3, 3, 3], + logical_constraint_list=[ + m.lim['mixer'].name, + m.lim['reactor'].name, + m.lim['centrifuge'].name, + ], + direction_norm=k, + time_limit=time_limit, + ) + end_time = time.time() + elapsed_time = end_time - start_time + + print(result) + new_results={'Solver': solver, 'k': str(k), 'Objective': pe.value(m.obj), 'Time': elapsed_time} + + writer.writerow([new_results['Solver'], new_results['k'], new_results['Objective'], new_results['Time']]) + + print(new_results) + + print("objective: ", pe.value(m.obj)) + for k in m.k: + for j in m.j: + print(f"Y[{k},{j}] = {pe.value(m.Y[k, j])}")