3. PID design pipeline¶
In this example, the pipeline design described in the paper Merigo et al., 2019 is fully reproduced. All the code detail to perform the simulation is available on the GitHub project here.
[ ]:
# Third party imports
import numpy as np
import pandas as pd
from bokeh.models import HoverTool
from bokeh.layouts import row, column
from bokeh.plotting import figure, show
from bokeh.io import output_notebook
# Local imports
from pid_sim_functions import simu, Patient_table
from python_anesthesia_simulator import metrics
output_notebook()
3.1. Load controller parameters¶
Load parameters from .csv file
[ ]:
phase = 'maintenance'
# phase = 'induction'
if phase == 'maintenance':
param_opti = pd.read_csv('optimal_parameters_PID_reject.csv')
else:
param_opti = pd.read_csv('optimal_parameters_PID.csv')
3.2. Perform table population simulation¶
Test all the drug ratio on all the patient in the table
[3]:
ts = 1
IAE_list = []
if phase == 'induction':
TT_list = []
ST10_list = []
elif phase == 'maintenance':
TTp_list = []
TTn_list = []
p1 = figure(width=600, height=300)
p2 = figure(width=600, height=300)
p3 = figure(width=600, height=300)
for ratio in range(2, 3):
print('ratio = ' + str(ratio))
Kp = float(param_opti.loc[param_opti['ratio'] == ratio, 'Kp'].iloc[0])
Ti = float(param_opti.loc[param_opti['ratio'] == ratio, 'Ti'].iloc[0])
Td = float(param_opti.loc[param_opti['ratio'] == ratio, 'Td'].iloc[0])
PID_param = [Kp, Ti, Td, ratio]
for i in range(1, 14):
Patient_info = Patient_table[i-1][1:]
IAE, dataframe, BIS_param = simu(Patient_info, phase, PID_param)
source = pd.DataFrame(data=dataframe['BIS'], columns=['BIS'])
source.insert(len(source.columns), "time", dataframe['Time']/60)
source.insert(len(source.columns), "Ce50_P", BIS_param[0])
source.insert(len(source.columns), "Ce50_R", BIS_param[1])
source.insert(len(source.columns), "gamma", BIS_param[2])
source.insert(len(source.columns), "beta", BIS_param[3])
source.insert(len(source.columns), "E0", BIS_param[4])
source.insert(len(source.columns), "Emax", BIS_param[5])
plot = p1.line(x='time', y='BIS', source=source)
tooltips = [('Ce50_P', "@Ce50_P"), ('Ce50_R', "@Ce50_R"),
('gamma', "@gamma"), ('beta', "@beta"),
('E0', "@E0"), ('Emax', "@Emax")]
p1.add_tools(HoverTool(renderers=[plot], tooltips=tooltips))
p2.line(dataframe['Time']/60,dataframe['MAP'], legend_label='MAP (mmgh)')
p2.line(dataframe['Time']/60, dataframe['CO']*10,
legend_label='CO (cL/min)', line_color="#f46d43")
p3.line(dataframe['Time']/60, dataframe['u_propo'],
line_color="#006d43", legend_label='propofol (mg/min)')
p3.line(dataframe['Time']/60, dataframe['u_remi'],
line_color="#f46d43", legend_label='remifentanil (ng/min)')
metric = metrics.compute_control_metrics(dataframe['Time'].to_list(), dataframe['BIS'].to_list(), phase=phase)
if phase == 'induction':
TT_list.append(metric['TT'].iloc[0])
ST10_list.append(metric['ST10'].iloc[0])
elif phase == 'maintenance':
TTp_list.append(metric['TTp'].iloc[0])
TTn_list.append(metric['TTn'].iloc[0])
IAE_list.append(IAE)
p1.title.text = 'BIS'
p3.title.text = 'Infusion rates'
p3.xaxis.axis_label = 'Time (min)'
grid = row(column(p3, p1, p2))
show(grid)
print(f"Mean IAE : {np.mean(IAE_list):.2f}")
if phase == 'induction':
print(f"Mean TT : {np.mean(TT_list):.2f} minutes")
print(f"Min TT : {np.min(TT_list):.2f} minutes")
print(f"Max TT : {np.max(TT_list):.2f} minutes")
print(f"Mean ST10 : {np.round(np.nanmean(ST10_list),2):.2f} minutes")
print(f"Min ST10 : {np.round(np.nanmin(ST10_list),2):.2f} minutes")
print(f"Max ST10 : {np.round(np.nanmax(ST10_list),2):.2f} minutes")
elif phase == 'maintenance':
print(f"Mean TTp : {np.mean(TTp_list):.2f} minutes")
print(f"Min TTp : {np.min(TTp_list):.2f} minutes")
print(f"Max TTp : {np.max(TTp_list):.2f} minutes")
print(f"Mean TTn : {np.mean(TTn_list):.2f} minutes")
print(f"Min TTn : {np.min(TTn_list):.2f} minutes")
print(f"Max TTn : {np.max(TTn_list):.2f} minutes")
ratio = 2
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
Mean IAE : 1500.09
Mean TTp : 0.21 minutes
Min TTp : 0.13 minutes
Max TTp : 0.35 minutes
Mean TTn : 1.14 minutes
Min TTn : 0.82 minutes
Max TTn : 1.65 minutes
3.3. Inter-Patient Variability Test¶
Create a wide population with uniform distribution for patient demographie data and normal distribution for PD model
[4]:
#Simulation parameter
ratio = 2
Number_of_patient = 50
# Controller parameters
Kp = float(param_opti.loc[param_opti['ratio'] == ratio, 'Kp'].iloc[0])
Ti = float(param_opti.loc[param_opti['ratio'] == ratio, 'Ti'].iloc[0])
Td = float(param_opti.loc[param_opti['ratio'] == ratio, 'Td'].iloc[0])
PID_param = [Kp, Ti, Td, ratio]
IAE_list = []
if phase == 'induction':
TT_list = []
ST10_list = []
elif phase == 'maintenance':
TTp_list = []
TTn_list = []
p1 = figure(width=600, height=300)
p2 = figure(width=600, height=300)
p3 = figure(width=600, height=300)
for i in range(Number_of_patient):
#Generate random patient information with uniform distribution
age = np.random.randint(low=18,high=70)
height = np.random.randint(low=150,high=190)
weight = np.random.randint(low=50,high=100)
sex = np.random.randint(low=0,high=1)
Patient_info = [age, height, weight, sex] + [None]*6
IAE, dataframe, BIS_param = simu(Patient_info, phase, PID_param)
source = pd.DataFrame(data=dataframe['BIS'], columns=['BIS'])
source.insert(len(source.columns), "time", dataframe['Time']/60)
source.insert(len(source.columns), "Ce50_P", BIS_param[0])
source.insert(len(source.columns), "Ce50_R", BIS_param[1])
source.insert(len(source.columns), "gamma", BIS_param[2])
source.insert(len(source.columns), "beta", BIS_param[3])
source.insert(len(source.columns), "E0", BIS_param[4])
source.insert(len(source.columns), "Emax", BIS_param[5])
plot = p1.line(x='time', y='BIS', source=source)
tooltips = [('Ce50_P', "@Ce50_P"), ('Ce50_R', "@Ce50_R"),
('gamma', "@gamma"), ('beta', "@beta"),
('E0', "@E0"), ('Emax', "@Emax")]
p1.add_tools(HoverTool(renderers=[plot], tooltips=tooltips))
p2.line(dataframe['Time']/60,dataframe['MAP'], legend_label='MAP (mmgh)')
p2.line(dataframe['Time']/60, dataframe['CO']*10,
legend_label='CO (cL/min)', line_color="#f46d43")
p3.line(dataframe['Time']/60, dataframe['u_propo'],
line_color="#006d43", legend_label='propofol (mg/min)')
p3.line(dataframe['Time']/60, dataframe['u_remi'],
line_color="#f46d43", legend_label='remifentanil (ng/min)')
if phase == 'induction':
metric = metrics.compute_control_metrics(dataframe['Time'], dataframe['BIS'], phase=phase)
TT_list.append(metric['TT'].iloc[0])
ST10_list.append(metric['ST10'].iloc[0])
elif phase == 'maintenance':
metric = metrics.compute_control_metrics(dataframe['Time'], dataframe['BIS'], phase=phase)
TTp_list.append(metric['TTp'].iloc[0])
TTn_list.append(metric['TTn'].iloc[0])
IAE_list.append(IAE)
p1.title.text = 'BIS'
p3.title.text = 'Infusion rates'
p3.xaxis.axis_label = 'Time (min)'
grid = row(column(p3, p1, p2))
show(grid)
print(f"Mean IAE : {np.mean(IAE_list):.2f}")
if phase == 'induction':
print(f"Mean TT : {np.mean(TT_list):.2f} minutes")
print(f"Min TT : {np.min(TT_list):.2f} minutes")
print(f"Max TT : {np.max(TT_list):.2f} minutes")
print(f"Mean ST10 : {np.round(np.nanmean(ST10_list),2):.2f} minutes")
print(f"Min ST10 : {np.round(np.nanmin(ST10_list),2):.2f} minutes")
print(f"Max ST10 : {np.round(np.nanmax(ST10_list),2):.2f} minutes")
elif phase == 'maintenance':
print(f"Mean TTp : {np.mean(TTp_list):.2f} minutes")
print(f"Min TTp : {np.min(TTp_list):.2f} minutes")
print(f"Max TTp : {np.max(TTp_list):.2f} minutes")
print(f"Mean TTn : {np.mean(TTn_list):.2f} minutes")
print(f"Min TTn : {np.min(TTn_list):.2f} minutes")
print(f"Max TTn : {np.max(TTn_list):.2f} minutes")