.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/advanced/SI_engine_optimization.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_advanced_SI_engine_optimization.py: .. _ref_SI_engine_optimization: =========================================================== Find the optimal operating parameters for a given SI engine =========================================================== Finding the optimal solution is one of the prominent parts of the research and development activities. When it comes to *PyChemkin*, reaction mechanism (rate parameters) optimization, fuel surrogate blend optimization, and process/reactor parameters optimization are some of the common applications. This example project describes the steps of integrating *PyChemkin* and *pymoo*, a well-known Python multi-objectives optimization package, to determine the optimal combination of the "start of combustion (SOC) timing" and the "exhaust gas recirculation (EGR) ratio" such that the engine would generate the highest power output with minimal nitric oxide (NO) emission and no engine knock. The engine power output in this example is measured by the indicated mean effective pressure (IMEP) between the intake valve close (IVC) and the exhaust valve open (EVO). The NO emission is indicated by the cylinder-averaged peak NO mass fraction during the same time period. The engine knock is defined as the occurrence of auto-ignition of the "endgas" which is the gas mixture of the "unburned zone" of the SI engine model. The intensity of the engine knock in this example is represented by the peak value of the total thermicity of the endgas. When the peak thermicity value of a particular engine run goes above the preset threshold, the occurrence of engine knock is predicted. Consequently, the operating condition and the solution associated with the run will be discarded because of the violation of the no engine knock constraint enforced by the optimization process. The multi-objective algorithm ``NGSA-II`` is employed to perform the optimization. The best pair of SOC and EGR ratio that satisfies the no engine knock constraint will be "determined" by applying a "decision-making" method from the *pymoo* package. The parameters and the solutions of cases around the Pareto frontier of the optimization process will be plotted in the "design space" and the "objective space", respectively, with the best parameter/solution marked. Details of the Chemkin 0-D SI engine model are available at the "Ansys Product Help" site. Information about *pymoo* can be found online at pymoo.org. .. GENERATED FROM PYTHON SOURCE LINES 61-63 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 65-67 Import PyChemkin packages and start the logger ============================================== .. GENERATED FROM PYTHON SOURCE LINES 67-107 .. code-block:: Python import datetime import multiprocessing import os from pathlib import Path import ansys.chemkin.core as ck # Chemkin # chemkin spark ignition (SI) engine model (transient) from ansys.chemkin.core.engines.SI import SIengine from ansys.chemkin.core.inlet import Stream # external gaseous inlet from ansys.chemkin.core.logger import logger from ansys.chemkin.core.utilities import WorkingFolders import matplotlib.pyplot as plt # plotting import numpy as np # number crunching # pymoo: multi-objective optimization in Python from pymoo.algorithms.moo.nsga2 import NSGA2 from pymoo.core.problem import ElementwiseProblem, StarmapParallelization from pymoo.decomposition.asf import ASF from pymoo.operators.sampling.lhs import LHS from pymoo.optimize import minimize # check working directory current_dir = str(Path.cwd()) # create a separate working folder for the optimization runs top_level = WorkingFolders("SI_optimize", current_dir) top_dir = top_level.work top_level.done() # logger.debug("working directory: " + str(Path.cwd())) # set verbose mode ck.set_verbose(False) # set interactive mode for plotting the results # interactive = True: display plot # interactive = False: save plot as a PNG file global interactive interactive = False .. GENERATED FROM PYTHON SOURCE LINES 108-113 Create an SI engine calculator class ===================================== Create a local class that wraps around the actual ``SIengine`` class to make the setup of the SI engine parameter study more convenient. .. GENERATED FROM PYTHON SOURCE LINES 113-301 .. code-block:: Python class SIengineCalculator: """SI engine calculator with fixed set up parameters.""" def __init__(self, fresh_mixture: Stream): """Create an SI engine calculator.""" """ Create an SI engine calculator that instantiates an SI engine object with the given fresh (unburnt) mixture condition and predefined engine parameters. Parameters ---------- fresh_mixture: Mixture object the initial/fresh/unburnt condition """ # instantiate an SI engine object # set up the run and working directory name self.name = "SI_engine" self.si_calculator = SIengine(fresh_mixture, self.name) # run status self.runstatus = -100 # IMEP (objective #1) self.imep = 0.0 # peak NO mass fraction (objective #2) self.max_no = 1.0 # calculated maximum thermicity value in the unburned zone # [1/sec] (objective #3) self.max_sigma = 1.0e5 # peak endgas temperature [K] self.max_temp = 0.0 # set the required premixed flame model parameters self.set_engine_parameters() # set burn profile (variables/parameters) # the burn profile will be set right before each run during # the optimization process def run(self) -> int: """Run the SI engine calculation.""" """ Run individual SI engine calculation. Returns ------- status: integer run status code """ # reset status self.runstatus = -100 # run the flame speed calculation self.runstatus = self.si_calculator.run() # extract the laminar flame speed from the solution if self.runstatus == 0: # postprocess the unburned zone solution unburnedzone = 1 # burnedzone = 2 # because the memory is shared, it must be done as soon as # the run is finished # get solution variables from the unburned zone ierr = self.si_calculator.process_engine_solution(zone_id=unburnedzone) if ierr == 0: temperature = self.si_calculator.get_solution_variable_profile( "temperature" ) thermicity = self.si_calculator.get_solution_variable_profile( "thermicity" ) self.imep = self.si_calculator.get_engine_imep() self.max_temp = np.max(temperature) self.max_sigma = np.max(thermicity) # get the cylinder-averaged solution variables ierr = self.si_calculator.process_average_engine_solution() no_mass_frac = self.si_calculator.get_solution_variable_profile("no") self.max_no = np.max(no_mass_frac) del thermicity, temperature, no_mass_frac self.runstatus = ierr else: self.runstatus = ierr return self.runstatus def set_engine_parameters(self): """Set up basic engine parameters.""" """ Set up basic engine parameters related to key geometries and operating conditions. """ # Set up basic engine parameters # These engine parameters are used to describe the cylinder volume during the # simulation. The ``starting_CA`` argument should be the crank angle # corresponding to the cylinder IVC. The ``ending_CA`` argument is typically # the EVO crank angle. cylinder bore diameter [cm] self.si_calculator.bore = 8.5 # engine stroke [cm] self.si_calculator.stroke = 10.82 # connecting rod length [cm] self.si_calculator.connecting_rod_length = 17.853 # compression ratio [-] self.si_calculator.compression_ratio = 12 # engine speed [RPM] self.si_calculator.rpm = 900 # simulation start CA [degree] self.si_calculator.starting_ca = -120.2 # simulation end CA [degree] self.si_calculator.ending_ca = 139.8 # wall heat transfer parameters # Set up engine wall heat transfer model heattransferparameters = [0.15, 0.8, 0.0] # set cylinder wall temperature [K] t_wall = 434.0 self.si_calculator.set_wall_heat_transfer( "dimensionless", heattransferparameters, t_wall ) # in-cylinder gas velocity correlation parameter (Woschni) # [ ] gv_parameters = [2.28, 0.318, 0.324, 0.0] self.si_calculator.set_gas_velocity_correlation(gv_parameters) # set piston head top surface area [cm2] self.si_calculator.set_piston_head_area(area=56.75) # set cylinder clearance surface area [cm2] self.si_calculator.set_cylinder_head_area(area=56.75) # other engine model parameters # set tolerances in tuple: (absolute tolerance, relative tolerance) self.si_calculator.tolerances = (1.0e-12, 1.0e-8) # turn on the force non-negative solutions option in the solver self.si_calculator.force_nonnegative = False # set minimum zonal mass [g] self.si_calculator.set_minimum_zone_mass(minmass=1.0e-5) # set the maximum solver time step self.si_calculator.max_ca_step = 0.1 # set the number of crank angles between saving solution self.si_calculator.ca_step_for_saving_solution = 0.5 # set the number of crank angles between printing solution self.si_calculator.ca_step_for_printing_solution = 10.0 # turn ON adaptive solution saving self.si_calculator.adaptive_solution_saving(mode=True, steps=20) # set to detect auto-ignition # self.si_calculator.set_ignition_delay(method="T_inflection") # set species lower bound for equilibrium calculation self.si_calculator.set_burned_products_minimum_mole_fraction(1.0e-12) # suppress text output from the SI engine model self.si_calculator.stop_output() def set_burn_profile( self, start_combustion: float, burn_duration: float, wiebe_n: float, wiebe_b: float, ): """Set up the fuel burn rate profile.""" """ Set up the fuel burn rate profile using the Wiebe function. Parameters ---------- start_combustion: double start of combustion (SOC) crank angle [CA] burn_duration: double burn duration in crank angle [CA] wiebe_n: double Wiebe function n parameter value wiebe_b: double Wiebe function b parameter value """ # start of combustion CA self.si_calculator.set_burn_timing( soc=start_combustion, duration=burn_duration, ) self.si_calculator.wiebe_parameters(n=wiebe_n, b=wiebe_b) def get_solution_parameters(self) -> tuple[float, float, float]: """Engine performance parameters.""" """ Get the SI engine performance parameters to be optimized from the solution. Returns ------- imep: double indicated mean effective pressure [bar] max_no: double peak cylinder-averaged NO mass fraction [-] max_sigma: double peak total thermicity of the eng gas [1/sec] """ return self.imep, self.max_no, self.max_sigma .. GENERATED FROM PYTHON SOURCE LINES 302-307 Create an instance of the Chemistry Set ======================================== For PRF, the encrypted 14-component gasoline mechanism, ``gasoline_14comp_WBencrypted.inp``, is used. The chemistry set is named ``gasoline``. .. GENERATED FROM PYTHON SOURCE LINES 307-448 .. code-block:: Python def setup_fresh_mixture( fuel_compo: list[tuple[str, float]], phi: float, egr_rate: float ) -> Stream: """Set up the chemistry set and create the unburned mixture.""" """ Set up the chemistry set and create the fresh fuel-air mixture at the intake valve closing (IVC). Parameters ---------- fuel: tuple (str, double) fuel species mole fractions phi: double equivalence ratio egr_rate: double EGR rate ratio Returns ------- fresh: Mixture object fresh unburned mixture (initial charge) """ # set mechanism directory (the default Chemkin mechanism data directory) data_dir = Path(ck.ansys_dir) / "reaction" / "data" mechanism_dir = data_dir # create a chemistry set based on the gasoline 14 components mechanism MyGasMech = ck.Chemistry(label="Gasoline") # set mechanism input files # including the full file path is recommended MyGasMech.chemfile = str(mechanism_dir / "gasoline_14comp_WBencrypt.inp") ####################################### # Preprocess the gasoline chemistry set # ===================================== # preprocess the mechanism files _ = MyGasMech.preprocess() ################################################ # Set up the stoichiometric gasoline-air mixture # ============================================== # You must set up the stoichiometric gasoline-air mixture for the subsequent # SI engine calculations. Here the ``x_by_equivalence_ratio()`` method is used. # You create the ``fuel`` and the ``air`` mixtures first. You then define the # *complete combustion product species* and provide the *additives* composition # if applicable. # Finally, you simply set ``equivalenceratio=1`` to create the stoichiometric # gasoline-air mixture. # # For PRF 90 gasoline, the recipe is ``[("ic8h18", 0.9), ("nc7h16", 0.1)]``. # # .. note:: # The "surrogate blend utility" from the Chemkin "reaction workbench" can be # used to create a surrogate that matches the properties (heating value, # distillation curve, density, viscosity, etc.) of an actual fuel. # # create the fuel mixture fuelmixture = ck.Mixture(MyGasMech) # set fuel composition fuelmixture.x = fuel_compo # setting pressure and temperature is not required in this case fuelmixture.pressure = 1.15 * ck.P_ATM # Update the fresh mixture temperature according to the EGR ratio # intake gas temperature [K] temp_intake = 310.0 # engine exhaust gas temperature [K] temp_egr = 950.0 # gas charge temperature [K] fuelmixture.temperature = (1.0 - egr_rate) * temp_intake + egr_rate * temp_egr # create the oxidizer mixture: air air = ck.Mixture(MyGasMech) air.x = [("o2", 0.21), ("n2", 0.79)] # setting pressure and temperature is not required in this case air.pressure = fuelmixture.pressure air.temperature = fuelmixture.temperature # products from the complete combustion of the fuel mixture and air products = ["co2", "h2o", "n2"] # species mole fractions of added/inert mixture. # You can also create an additives mixture here. add_frac = np.zeros(MyGasMech.kk, dtype=np.double) # no additives: all zeros # create the unburned fuel-air mixture fresh = ck.Mixture(MyGasMech) # mean equivalence ratio equiv = phi iError = fresh.x_by_equivalence_ratio( MyGasMech, fuelmixture.x, air.x, add_frac, products, equivalenceratio=equiv ) ########################################################## # Specify pressure and temperature of the fuel-air mixture # ======================================================== # Since you are going to use ``fresh`` fuel-air mixture to instantiate # the engine object later, setting the mixture pressure and temperature # is equivalent to setting the initial temperature and pressure of the # engine cylinder. fresh.pressure = fuelmixture.pressure fresh.temperature = fuelmixture.temperature ########################################### # Add EGR to the fresh fuel-air mixture # ========================================= # Many engines have the configuration for exhaust gas recirculation (EGR). Chemkin # engine models let you add the EGR mixture to the fresh fuel-air mixture entering # the cylinder. If the engine you are modeling has EGR, you should have # the EGR ratio, which is generally the volume ratio of the EGR mixture and # the fresh fuel-air mixture. # However, because you know nothing about the composition of the exhaust gas, # you cannot simply combine these two mixtures. In this case, you use the # ``get_egr_mole_fraction()`` method to estimate the major components of # the exhaust gas from the combustion of the fresh fuel-air mixture. The # ``threshold=1.0e-8`` parameter tells the method to ignore any species with # a mole fraction below the threshold value. Once you have the EGR mixture # composition, use the ``x_by_equivalence_ratio()`` method a second time to # re-create the fuel-air mixture ``fresh`` with the original ``fuelmixture`` and # ``air`` mixtures, along with the EGR composition that you just got as the # *additives*. egr_ratio = egr_rate # compute the EGR stream composition in mole fractions add_frac = fresh.get_egr_mole_fraction(egr_ratio, threshold=1.0e-8) # recreate the initial mixture with EGR ierror = fresh.x_by_equivalence_ratio( MyGasMech, fuelmixture.x, air.x, add_frac, products, equivalenceratio=equiv, threshold=1.0e-8, ) if ierror != 0: print("error...creating the initial fuel-oxidizer mixture.") return fresh .. GENERATED FROM PYTHON SOURCE LINES 449-475 Set up the optimization problem ================================ Define the optimization problem. Set up the optimization parameters, variables, objectives, constraints, and the bounds of the variables. Since the ``_evaluate()`` method perform one SI engine simulation at a time, the optimization problem is defined as an ``ElementwiseProblem``. By implementing the SI engine model as an element-wise problem, the parallelization available in *pymoo* can be directly utilized. There are two design parameters/variables in this project. The variable SOC is used to set up the mass burned fraction profile of the SI engine model. And the EGR ratio is applied to create the initial gas charge inside the engine cylinder. The objective and the constraint values come from the result of the SI engine simulation. The first objective is the negative value of IMEP, the second objective is the peak NO mass fraction in the cylinder, and the constraint value is the knock intensity as indicated by the total thermicity of the endgas. In *pymoo*, the objectives are "minimized" so the objective of maximizing the IMEP must be reformulated into minimizing the negative value of IMEP, and the constraint must be written in the form of "g <= 0.0". .. note:: See the *pymoo* documentation for details about how the optimization problem is set up and evaluated. .. GENERATED FROM PYTHON SOURCE LINES 475-610 .. code-block:: Python class SI_engine_opt(ElementwiseProblem): """SI engine optimization object.""" def __init__( self, numb_vars: int, numb_objs: int, fuel: list[tuple[str, float]], phi: float, **kwargs, ): """Set up the SI engine optimization.""" """ Set up the SI engine optimization. The objective is to find the optimal fuel mass burn profile to maximize the engine indicated mean effective pressure (IMEP) while avoiding engine knock and minimizing the NO emission. The fuel mass burn profile is described by the anchor points: crank angles at 10%, 50%, and 90% mass burned. Parameters ---------- numb_vars: integer number of variables numb_objs: integer number of objectives fuel: tuple (str, double) fuel species mole fractions phi: double equivalence ratio of the initial gas charge """ # number of inequality constraints numb_constrs = 1 # lower bounds for the variables lower_bounds = np.zeros(numb_vars, dtype=np.double) upper_bounds = np.zeros(numb_vars, dtype=np.double) lower_bounds[0] = -15.0 lower_bounds[1] = 0.0 # upper bounds for the variables upper_bounds[0] = 10.0 upper_bounds[1] = 0.60 # initialize the "Problem" object super().__init__( n_var=numb_vars, n_obj=numb_objs, n_ieq_constr=numb_constrs, n_eq_constr=0, xl=lower_bounds, xu=upper_bounds, ) # number of evaluation calls self.call_count = 0 self.fail_count = 0 self.failed_cases: list[int] = [] # initial gas charge parameters # fuel composition self.fuel = fuel # equivalence ratio self.phi = phi def _evaluate(self, x, out, *args, **kwargs): """Evaluate the objective function values.""" """ Evaluate the objective function values with the given set of variable values. Parameters ---------- x: list[float] values of the pymoo variables to be optimized. out: dict[str: list[float]] pymoo OUT dictionary containing evaluated function values and constrain values """ # initialize "bad" return values to drive the search away from this state point f1 = 0.0 f2 = 1.0 # g1 = 5.0e4 # variables # x[0]: start of combustion crank angle # x[1]: exhaust gas recirculation (EGR) rate # create and change to the working directory for this run name = "SI_engine_" + str(self.call_count) sub_work_folder = WorkingFolders(name, top_dir) # set up the unburned mixture of the SI engine at IVC init_mixture = setup_fresh_mixture(self.fuel, self.phi, x[1]) # instantiate the SI engine engine_case = SIengineCalculator(init_mixture) # set the burned mass fraction profile # use two sets of Wiebe function parameters to describe # the burned mass fraction profile in the SI engine # profile #1: fast burn wiebe_n_egr = 1.8 wiebe_b_egr = 8.5 # profile #2: regular burn wiebe_n_fresh = 2.0 wiebe_b_fresh = 8.0 # EGR ratio r = x[1] r1 = 1.0 - r # burned mass profile as a function of EGR ratio wiebe_b = r * wiebe_b_egr + r1 * wiebe_b_fresh wiebe_n = r * wiebe_n_egr + r1 * wiebe_n_fresh burn_duration = 45.5 # degrees # set up the mass burned profile in the SI engine # the burn duration is fixed engine_case.set_burn_profile(x[0], burn_duration, wiebe_n, wiebe_b) # objects # f1: IMEP [bar] (to be maximized) # f2: peak cylinder-averaged NO mass fraction [-] (to be minimized) runstatus = engine_case.run() if runstatus == 0: # get the results only when the simulation is successful f1, f2, max_sigma = engine_case.get_solution_parameters() # compute the constraints # g1: minimal engine knock # max_sigma: peak endgas total thermicity (knock indicator) [1/sec] g1 = max_sigma - 5.0e4 else: # tally the number of failed runs self.fail_count += 1 self.failed_cases.append(self.call_count) # clean up del engine_case sub_work_folder.done() # self.call_count += 1 # return solution values to the optimizer # return negative solution value for maximization objective out["F"] = [-f1, f2] # return constraints out["G"] = [g1] .. GENERATED FROM PYTHON SOURCE LINES 611-621 Run the optimization ==================== Select the optimization algorithm and the decision making method. Since this project has two objectives and one constraint, the multi-objective algorithm ``NSGA-II`` is selected as the optimization algorithm. A relatively small population size of 50 is used for this simple optimization problem. The ``StarmapParallelization()`` method from *pymoo* is used to parallelize the SI engine simulation runs during the optimization process. .. GENERATED FROM PYTHON SOURCE LINES 621-791 .. code-block:: Python def run_optimization( numb_vars: int, numb_objs: int, fuel: list[tuple[str, float]], phi: float ): """Run the optimization process.""" """ Run the SI engine optimization process with pymoo. Parameters ---------- numb_vars: integer number of variables numb_objs: integer number of objectives fuel: tuple (str, double) fuel species mole fractions phi: double equivalence ratio of the initial gas charge Returns ------- res: pymoo Result object pymoo minimization results and history """ # create the "problem" to be optimized # set up the optimization algorithm # NSGA: Non-dominated Sorting Genetic Algorithm # LHS: Latin Hypercube Sampling algorithm = NSGA2( pop_size=50, sampling=LHS(), eliminate_duplications=True, ) # set up the thermination condition # run the optimization run_parallel = True if run_parallel: # parallel mode # initialize the process pool and create the runner # number of available cpu cores numb_cores = max(os.cpu_count(), 1) numb_workers = numb_cores - 2 numb_workers = 4 # max(numb_workers, 1) pool = multiprocessing.Pool(numb_workers) runner = StarmapParallelization(pool.starmap) # define the problem by passing the starmap interface of the thread pool problem = SI_engine_opt( numb_vars=numb_vars, numb_objs=numb_objs, fuel=fuel, phi=phi, elementwise_runner=runner, ) else: # serial mode problem = SI_engine_opt( numb_vars=numb_vars, numb_objs=numb_objs, fuel=fuel, phi=phi ) res = minimize(problem, algorithm, termination=("n_gen", 1), seed=1) # run summary # execution time [sec] print(f"Total wall time: {datetime.timedelta(seconds=res.exec_time)}") # run statistics print(f"Failed runs: {problem.fail_count}/{problem.call_count}") if problem.fail_count > 0: print("failed cases") for n in problem.failed_cases: print(f"run # {n}") ############################ # Post-process the solutions # ========================== # Parse the solutions from the optimization process and display them with # the best solution marked by a green cross mark. The variable values are stored # in the ``Result.X`` array and the corresponding objective values are in # the ``Result.F`` array. Because the objective of maximizing the "IMEP" has to be # converted to a minimization objective in *pymoo*, an alternative objective # of minimizing "-IMEP" is adopted. Consequently, the raw -IMEP objective solution # must be converted back to IMEP to make sense of the optimization and the # subsequent "decision-making" results. # # In this project, the decomposition function ``ASF()`` is applied to retrieve # the best solution. Weight factors are used to make the IMEP objective slightly # more important than the NO emission objective. # process the results x = res.X f = res.F # convert the objective result back to the original IMEP imep = -1.0 * f[:, 0] # approx_ideal = f.min(axis=0) approx_nadir = f.max(axis=0) # finding the optimal result # normalize the results diff = approx_ideal - approx_nadir if np.isclose(diff[0], 0.0, atol=1.0e-9) or np.isclose(diff[1], 0.0, atol=1.0e-9): nf = f - approx_ideal else: nf = (f - approx_ideal) / (approx_ideal - approx_nadir) # set the weight factors weight = np.zeros(n_objs, dtype=np.double) weight[0] = 1.0 / 0.6 weight[1] = 1.0 / 0.4 # find the best solution # use the Augmented Scalarization Function (ASF) method decomp = ASF() i = decomp.do(nf, weight).argmin() print("Best outcome according to ASF:") print(f"run # = {i}") print(f"Predicted IMEP = {imep[i]} [bar]") print(f"Predicted peak NO mass fraction = {f[i, 1]} [-]") print(f"Start of combustion = {x[i, 0]} [degree CA ATDC]") print(f"EGR ratio = {x[i, 1] * 100.0} % vol.") # plot the design space xl, xu = problem.bounds() plt.figure(figsize=(7, 5)) plt.title("Design Space") plt.scatter(x[:, 0], x[:, 1], s=30, facecolors="none", edgecolors="b") plt.scatter(x[i, 0], x[i, 1], color="green", marker="x", s=75, label="Best") plt.xlim(xl[0], xu[0]) plt.ylim(xl[1], xu[1]) plt.ylabel("EGR ratio [-]") plt.xlabel("Start of Combustion [degree CA]") plt.legend() if interactive: plt.show() else: plt.savefig("plot_SI_engine_optimization_designspace.png", bbox_inches="tight") # plot the objective space to show the Pareto frontier plt.figure(figsize=(7, 5)) plt.scatter(imep, f[:, 1], s=30, facecolors="none", edgecolors="b") # note that because IMEP = -F[:, 0], approx_ideal[0] and approx_nadir[0] # must be modified accordingly ideal_imep = -1.0 * approx_ideal[0] nadir_imep = -1.0 * approx_nadir[0] plt.scatter( ideal_imep, approx_ideal[1], facecolors="none", edgecolors="red", marker="*", s=100, label="Ideal Point (Approx)", ) plt.scatter( nadir_imep, approx_nadir[1], facecolors="none", edgecolors="black", marker="p", s=100, label="Nadir Point (Approx)", ) plt.scatter(imep[i], f[i, 1], color="green", marker="x", s=75, label="Best") plt.title("Objective Space - Power & Emission") plt.xlabel("IMEP [bar]") plt.ylabel("Peak NO mass fraction [-]") plt.legend() if interactive: plt.show() else: plt.savefig("plot_SI_engine_optimization_objspace.png", bbox_inches="tight") # del imep, nf, xl, xu # return res .. GENERATED FROM PYTHON SOURCE LINES 792-802 Set up the main driver of the optimization process ================================================== A main method for the optimization process is necessary especially with the intention to run the design point cases in parallel. You can further post-process the results from the optimization process. Here the correlations between the design variables (SOC and EGR ratio) and the objectives (IMEP and NO emission) are explored. You can consult *pymoo* online API document to find out additional information offered by the ``Result`` object. .. GENERATED FROM PYTHON SOURCE LINES 802-845 .. code-block:: Python if __name__ == "__main__": # properties of the initial gas charge # fuel composition # the "surrogate blend utility" from the Chemkin "reaction workbench" can be # used to create a surrogate that matches the properties (heating value, # distillation curve, density, viscocity, etc.) of an actual fuel fuel = [("ic8h18", 0.6), ("nc7h16", 0.0), ("mch", 0.1), ("c6h5c3h7", 0.3)] # equivalence ratio phi = 1.0 # perform optimization # number of variables n_vars = 2 # number of objectives n_objs = 2 opt_results = run_optimization(n_vars, n_objs, fuel, phi) # process the optimization results x = opt_results.X f = opt_results.F imep = -1.0 * f[:, 0] # plot correlations between the variable/parameter and the objectives plt.figure(figsize=(7, 5)) plt.scatter(x[:, 0], imep, s=30, facecolors="none", edgecolors="r") plt.title("Correlation - SOC & Power") plt.ylabel("IMEP [bar]") plt.xlabel("Start of combustion [CA]") if interactive: plt.show() else: plt.savefig("plot_SI_engine_optimization_correlation1.png", bbox_inches="tight") # plt.figure(figsize=(7, 5)) plt.scatter(x[:, 1], imep, s=30, facecolors="none", edgecolors="r") plt.title("Correlation - EGR & Power") plt.ylabel("IMEP [bar]") plt.xlabel("EGR ratio [-]") # clean up ck.done() if interactive: plt.show() else: plt.savefig("plot_SI_engine_optimization_correlation2.png", bbox_inches="tight") .. _sphx_glr_download_examples_advanced_SI_engine_optimization.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: SI_engine_optimization.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: SI_engine_optimization.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: SI_engine_optimization.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_