from fmpy import * from fmpy import read_model_description, extract from fmpy.fmi2 import FMU2Slave import numpy as np import shutil import pandas as pd import random import plotly.graph_objects as go def add_noise(min, max, std, mean): return np.clip(np.random.normal(mean, std), min, max) def gen_pulsive_noise(count, min_val, max_val): pulse = 0 n = count // 100 # Use integer division to get the floor of count / 100 start = 100 * n end = start + 50 if (count > start and count < end and n % 15 == 0): pulse = random.randint(min_val, max_val) return pulse def simulation(input_noise, output_noise, min, max, std, mean, min_pul, max_pul, kp, ki, kd, show_original): vrs = {} fmu = 'Pharmacokinetics_4_comportmental_model_PI_ref_FMU_base4_OAAS_lnx.fmu' model_description = read_model_description(fmu) for variable in model_description.modelVariables: vrs[variable.name] = variable.valueReference start_time = 0.0 stop_time = 7000 step_size = 1 unzipdir = extract(fmu) fmu = FMU2Slave(guid=model_description.guid, unzipDirectory=unzipdir, modelIdentifier=model_description.coSimulation.modelIdentifier, instanceName='instance1') # initialize fmu.instantiate() fmu.setupExperiment(startTime=start_time) fmu.enterInitializationMode() fmu.exitInitializationMode() fmu.setReal([vrs["amesim_interface.Age_year"]], [60]) fmu.setReal([vrs["amesim_interface.BIS0"]], [95.6]) fmu.setReal([vrs["amesim_interface.BISmin"]], [8.9]) fmu.setReal([vrs["amesim_interface.Drug_concentration_mgmL"]], [20]) fmu.setReal([vrs["amesim_interface.EC50"]], [2.23]) fmu.setReal([vrs["amesim_interface.Gamma"]], [1.58]) fmu.setReal([vrs["amesim_interface.Gender_0male_1female"]], [1]) fmu.setReal([vrs["amesim_interface.Height_cm"]], [168]) fmu.setReal([vrs["amesim_interface.Infusion_rate_mLh"]], [200]) fmu.setReal([vrs["amesim_interface.Weight_kg"]], [75]) vr_input = vrs["amesim_interface.Infusion_rate_mLh"] vr_output = vrs["amesim_interface.BIS_Index"] rows = [] # list to record the results time = start_time infusion_rate = 200 i = 0 target = 30 last_error = 0 # simulation loop while time < stop_time: if time >= 2.4e3 and time < 4.5e3: target = 70 p = 0 i = 0 if time >= 4.5e3: target = 30 p = 0 i = 0 bis = fmu.getReal([int(vr_output)])[0] if time > step_size else 95.6 if input_noise: bis += add_noise(min, max, std, mean) p = bis - target i = i + p d = p - last_error last_error = p infusion_rate = np.clip(kp*p + ki*i + kd*d, 0, 200) if output_noise: infusion_rate += gen_pulsive_noise(time, min_pul, max_pul) fmu.setReal([vr_input], [int(infusion_rate)]) # perform one step fmu.doStep(currentCommunicationPoint=time, communicationStepSize=step_size) # advance the time time += step_size # get the values for 'inputs' and 'outputs[4]' inputs, outputs = fmu.getReal([int(vr_input), int(vr_output)]) # append the results rows.append((time, bis, inputs)) fmu.terminate() fmu.freeInstance() shutil.rmtree(unzipdir, ignore_errors=True) result = np.array(rows, dtype=np.dtype([('time', np.float64), ('BIS', np.float64), ('Infusion', np.float64)])) result_baseline = np.load("result.npy") df = pd.DataFrame(result) df_original = pd.DataFrame(result_baseline) trace1 = go.Scatter(x=df.index, y=df['BIS'], mode='lines', name='BIS') fig1 = go.Figure(data=trace1) fig1.update_layout(height=400, width=1200, title_text="BIS evolution") # Add a line trace for column_2 in the second subplot trace2 = go.Scatter(x=df.index, y=df['Infusion'], mode='lines', name='Infusion') fig2 = go.Figure(data=trace2) fig2.update_layout(height=400, width=1200, title_text="Infusion evolution") if show_original: fig1.add_trace(go.Scatter(x=df_original.index, y=df_original['BIS'], mode='lines', name='BIS original', line=dict(color="red"), opacity=0.5)) fig2.add_trace(go.Scatter(x=df_original.index, y=df_original['Infusion'], mode='lines', name='Infusion original', line=dict(color="red"), opacity=0.5)) return fig1, fig2 import gradio as gr with gr.Blocks() as demo: with gr.Row(): with gr.Column(scale=1): gr.Markdown("**BIS noise**") input_noise = gr.inputs.Checkbox(label="Add Gaussian noise") with gr.Accordion("noise range"): min_gaussian = gr.inputs.Slider(minimum=0, maximum=20, step=1, default=5, label="noise min") max_gaussian = gr.inputs.Slider(minimum=0, maximum=20, step=1, default=10, label="noise max") std_gaussian = gr.inputs.Slider(minimum=0, maximum=10, step=1, default=2, label="noise standard deviation") mean_gaussian = gr.inputs.Slider(minimum=0, maximum=20, step=1, default=7, label="noise mean") gr.Markdown("**Infusion noise**") with gr.Blocks(): output_noise = gr.inputs.Checkbox(label="Add Pulsive noise") with gr.Accordion("noise range"): min_pul = gr.inputs.Slider(minimum=0, maximum=50, step=1, default=50, label="noise min") max_pul = gr.inputs.Slider(minimum=0, maximum=150, step=1, default=150, label="noise max") gr.Markdown("**PID controller paramters**") with gr.Blocks(): kp_slider = gr.inputs.Slider(minimum=0, maximum=20, default=4, label="kp") ki_slider = gr.inputs.Slider(minimum=0, maximum=1, default=0.01, label="ki") kd_slider = gr.inputs.Slider(minimum=0, maximum=200, default=120, label="kd") show_original = gr.inputs.Checkbox(label="Show original") button = gr.Button("Simulate") with gr.Column(scale=5): plot1 = gr.Plot(label="BIS evolution") plot2 = gr.Plot(label="Infusion evolution") button.click(simulation, inputs=[input_noise, output_noise, min_gaussian, max_gaussian, std_gaussian, mean_gaussian, min_pul, max_pul, kp_slider, ki_slider, kd_slider, show_original], outputs=[plot1, plot2]) demo.launch()