From dac3e3b61bd0b64749d1fd669c3c6ddcf12b5787 Mon Sep 17 00:00:00 2001 From: Joseph Hellerstein Date: Thu, 16 May 2024 14:00:47 -0700 Subject: [PATCH] Enable running multiple simulations of the same design parameters to address random noise; run all points in parallell. --- src/controlSBML/control_sbml.py | 4 ++ src/controlSBML/sbml_system.py | 4 +- src/controlSBML/siso_closed_loop_designer.py | 60 +++++++++++++------- tests/test_control_sbml.py | 2 +- 4 files changed, 46 insertions(+), 24 deletions(-) diff --git a/src/controlSBML/control_sbml.py b/src/controlSBML/control_sbml.py index 636b802..268d1c1 100644 --- a/src/controlSBML/control_sbml.py +++ b/src/controlSBML/control_sbml.py @@ -653,6 +653,10 @@ def _plotClosedLoop(self, kP=kP, kI=kI, kD=kD, kF=kF, setpoint=setpoint, selections=selections, times=times, ) + if response_ts is None: + msg = "System is unstable for kP={}, kI={}, kD={}, kF={}".format(kP, kI, kD, kF) + msgs.warn(msg) + return response_ts, builder self._sbml_system.plotSISOClosedLoop(response_ts, setpoint, times=times, **plot_dct) return response_ts, builder diff --git a/src/controlSBML/sbml_system.py b/src/controlSBML/sbml_system.py index 20d3854..0504899 100644 --- a/src/controlSBML/sbml_system.py +++ b/src/controlSBML/sbml_system.py @@ -516,14 +516,13 @@ def getName(self, is_input=True): return names[0] def plotSISOClosedLoop(self, timeseries:Timeseries, setpoint, - selections=None, mgr:Optional[OptionManager]=None, **kwargs): + mgr:Optional[OptionManager]=None, **kwargs): """ Plots the results of a closed lop simulation. Input and output are defined in the SBMLSystem constructor. Args: timeseries: Timeseries setpoint: float - selections: list-str (names of species to be selected for plot. output_name is included by default.) kwargs: dict (kwargs for plotOneTS) """ input_name = self.getName(is_input=True) @@ -532,7 +531,6 @@ def plotSISOClosedLoop(self, timeseries:Timeseries, setpoint, if mgr is None: mgr = OptionManager(new_kwargs) new_kwargs["is_plot"] = False - import pdb; pdb.set_trace() df = pd.DataFrame(timeseries[output_name], columns=[output_name]) new_kwargs.setdefault(cn.O_AX2, 0) plot_result = util.plotOneTS(df, colors=[cn.SIMULATED_COLOR], **new_kwargs) diff --git a/src/controlSBML/siso_closed_loop_designer.py b/src/controlSBML/siso_closed_loop_designer.py index f730a04..2b4b98c 100644 --- a/src/controlSBML/siso_closed_loop_designer.py +++ b/src/controlSBML/siso_closed_loop_designer.py @@ -38,7 +38,10 @@ Workunit = collections.namedtuple("Workunit", "system input_name output_name setpoint times is_greedy num_restart is_report") -DesignResult = collections.namedtuple("DesignResult", "dataframe") +# DesignResult +# dataframe: table of design results +# max_count: maximum count of the design parameters that successfully simulate +DesignResult = collections.namedtuple("DesignResult", "dataframe max_count") def _calculateClosedLoopTransferFunction(open_loop_transfer_function=None, kP=None, kI=None, kD=None, kF=None, @@ -291,6 +294,7 @@ def designAlongGrid(self, grid:Grid, is_greedy:bool=False, num_restart:int=1, is_report:bool=False, num_process:int=-1)->DesignResult: """ Design objective: Create a feasible system (stable, no negative inputs/outputs) that minimizes residuals. + For systems with random noise or disturbance, uses the maximum value of the score. Args: grid: Grid @@ -309,34 +313,50 @@ def getValue(val): # point_evaluator = PointEvaluator(self.system.copy(), self.input_name, self.output_name, self.setpoint, self.sign, self.times, is_greedy=is_greedy) - parallel_search = ParallelSearch(point_evaluator, grid.points, num_process=num_process, is_report=is_report) # type: ignore - search_results = [] + # Expand the number of points + points = [] for _ in range(num_restart): - parallel_search.search() - search_results.append(parallel_search.getSearchResults()) - if len(search_results) == 0: + points.extend(grid.points) + # Do the search + parallel_search = ParallelSearch(point_evaluator, points, num_process=num_process, is_report=is_report) # type: ignore + parallel_search.search() + search_result_df = parallel_search.getSearchResults() + if len(search_result_df) == 0: df = pd.DataFrame([[None, None, None, None, cn.DESIGN_RESULT_CANNOT_SIMULATE]], columns=[CP_kP, CP_kI, CP_kD, CP_kF, cn.SCORE, cn.REASON]) - return DesignResult(dataframe=df) + return DesignResult(dataframe=df, max_count=0) # Merge the results and sort by score - search_result_df = pd.concat(search_results) search_result_df = search_result_df.sort_values(cn.SCORE, ascending=True) # type: ignore search_result_df = search_result_df.reset_index(drop=True) + # Handle replications of the same design parameters and select successful designs + search_result_df = search_result_df[search_result_df[cn.SCORE].notna()] + del search_result_df[cn.REASON] + sort_columns = list(grid.axis_dct.keys()) + groupby = search_result_df.groupby(sort_columns) + count_df = groupby.count() + max_count = np.max(count_df[cn.SCORE]) + sel = count_df == max_count + mean_df = groupby.max() + threshold_mean_df = mean_df[sel] + threshold_mean_df = threshold_mean_df[threshold_mean_df[cn.SCORE].notna()] + sorted_mean_df = threshold_mean_df.reset_index() + sorted_mean_df = sorted_mean_df.sort_values(cn.SCORE, ascending=True) + sorted_mean_df = sorted_mean_df.reset_index() # Record the result - self.residual_mse = search_result_df.loc[0, cn.SCORE] - if CP_kP in search_result_df.columns: - self.kP = getValue(search_result_df.loc[0, CP_kP]) - if CP_kI in search_result_df.columns: - self.kI = getValue(search_result_df.loc[0, CP_kI]) - if CP_kD in search_result_df.columns: - self.kD = getValue(search_result_df.loc[0, CP_kD]) - if CP_kF in search_result_df.columns: - self.kF = getValue(search_result_df.loc[0, CP_kF]) + self.residual_mse = sorted_mean_df.loc[0, cn.SCORE] + if CP_kP in sorted_mean_df.columns: + self.kP = getValue(sorted_mean_df.loc[0, CP_kP]) + if CP_kI in sorted_mean_df.columns: + self.kI = getValue(sorted_mean_df.loc[0, CP_kI]) + if CP_kD in sorted_mean_df.columns: + self.kD = getValue(sorted_mean_df.loc[0, CP_kD]) + if CP_kF in sorted_mean_df.columns: + self.kF = getValue(sorted_mean_df.loc[0, CP_kF]) # Save the results if self.save_path is not None: - search_result_df.to_csv(self.save_path, index=False) - self._design_result_df = search_result_df.copy() - return DesignResult(dataframe=search_result_df) + sorted_mean_df.to_csv(self.save_path, index=False) + self._design_result_df = sorted_mean_df + return DesignResult(dataframe=sorted_mean_df, max_count=max_count) @property def design_result_df(self): diff --git a/tests/test_control_sbml.py b/tests/test_control_sbml.py index d3b33e3..ba33169 100644 --- a/tests/test_control_sbml.py +++ b/tests/test_control_sbml.py @@ -497,7 +497,7 @@ def testBug15(self): grid.addAxis("kP", min_value=0.5, max_value=10, num_coordinate=10) grid.addAxis("kI", min_value=0.002, max_value=0.02, num_coordinate=10) SETPOINT = 1000 - design_result = ctlsb.plotGridDesign(grid, times=TIMES, setpoint=SETPOINT) + design_result = ctlsb.plotGridDesign(grid, times=TIMES, setpoint=SETPOINT, num_restart=10)