3. PID design pipeline
In this example, the pipeline design described in the paper [Merigo2019]_ is fully reproduced. All the code detail to perform the simulation is available on the github project here
[Merigo2019]
Merigo, F. Padula, N. Latronico, M. Paltenghi, and A. Visioli, “Optimized PID control of propofol and remifentanil coadministration for general anesthesia,” Communications in Nonlinear Science and Numerical Simulation, vol. 72, pp. 194–212, Jun. 2019, doi: 10.1016/j.cnsns.2018.12.015.
[1]:
# 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
[2]:
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'])
Ti = float(param_opti.loc[param_opti['ratio'] == ratio, 'Ti'])
Td = float(param_opti.loc[param_opti['ratio'] == ratio, 'Td'])
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)')
if phase == 'induction':
TT, BIS_NADIR, ST10, ST20, US = metrics.compute_control_metrics(dataframe['Time'], dataframe['BIS'], phase=phase)
TT_list.append(TT)
ST10_list.append(ST10)
elif phase == 'maintenance':
TTp, BIS_NADIRp, TTn, BIS_NADIRn = metrics.compute_control_metrics(dataframe['Time'], dataframe['BIS'], phase=phase)
TTp_list.append(TTp)
TTn_list.append(TTn)
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("Mean IAE : " + str(np.mean(IAE_list)))
if phase == 'induction':
print("Mean TT : " + str(np.mean(TT_list)))
print("Min TT : " + str(np.min(TT_list)))
print("Max TT : " + str(np.max(TT_list)))
print("Mean ST10 : " + str(np.round(np.nanmean(ST10_list),2)))
print("Min ST10 : " + str(np.round(np.nanmin(ST10_list),2)))
print("Max ST10 : " + str(np.round(np.nanmax(ST10_list),2)))
elif phase == 'maintenance':
print("Mean TTp : " + str(np.mean(TTp_list)))
print("Min TTp : " + str(np.min(TTp_list)))
print("Max TTp : " + str(np.max(TTp_list)))
print("Mean TTn : " + str(np.mean(TTn_list)))
print("Min TTn : " + str(np.min(TTn_list)))
print("Max TTn : " + str(np.max(TTn_list)))
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.0851518615077
Mean TTp : 0.0
Min TTp : 0.0
Max TTp : 0.0
Mean TTn : 0.0
Min TTn : 0.0
Max TTn : 0.0
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'])
Ti = float(param_opti.loc[param_opti['ratio']==ratio, 'Ti'])
Td = float(param_opti.loc[param_opti['ratio']==ratio, 'Td'])
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)
gender = np.random.randint(low=0,high=1)
Patient_info = [age, height, weight, gender] + [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':
TT, BIS_NADIR, ST10, ST20, US = metrics.compute_control_metrics(dataframe['Time'], dataframe['BIS'], phase=phase)
TT_list.append(TT)
ST10_list.append(ST10)
elif phase == 'maintenance':
TTp, BIS_NADIRp, TTn, BIS_NADIRn = metrics.compute_control_metrics(dataframe['Time'], dataframe['BIS'], phase=phase)
TTp_list.append(TTp)
TTn_list.append(TTn)
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("Mean IAE : " + str(np.mean(IAE_list)))
if phase == 'induction':
print("Mean TT : " + str(np.mean(TT_list)))
print("Min TT : " + str(np.min(TT_list)))
print("Max TT : " + str(np.max(TT_list)))
print("Mean ST10 : " + str(np.round(np.nanmean(ST10_list),2)))
print("Min ST10 : " + str(np.round(np.nanmin(ST10_list),2)))
print("Max ST10 : " + str(np.round(np.nanmax(ST10_list),2)))
elif phase == 'maintenance':
print("Mean TTp : " + str(np.mean(TTp_list)))
print("Min TTp : " + str(np.min(TTp_list)))
print("Max TTp : " + str(np.max(TTp_list)))
print("Mean TTn : " + str(np.mean(TTn_list)))
print("Min TTn : " + str(np.min(TTn_list)))
print("Max TTn : " + str(np.max(TTn_list)))