Please check and edit my Python code.

$$Create a graphical output of a mathematical model using arbitrary variables and Python code, but you must use a PID controller.

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import odeint

# Constants

R $=287\mathrm{.}1$ # Gas constant $[$Nm$/$kgK$]$

T $=293$ # Temperature inside the chamber $[$K$]$

V $=0\mathrm{.}8$ # Volume of the chamber $[$L$]$

k $=1\mathrm{.}4$ # Adiabatic coefficient

p$\_$atm $=101325$ # Atmospheric pressure $[$Pa$]$

S$\_$v $=15\mathrm{.}1$ # Nominal pumping speed of the vacuum pump $[$m$^3/$h$]$

p$\_$v $=0\mathrm{.}7$ # Ultimate gas pressure of the vacuum pump $[$Pa$]$

V$\_$t $=1\mathrm{.}3$ # Volume inside the bellows $[$L$]$

# Variables

D $=1$ # Diameter of valve port

A$\_$cu $=1$ # Effective cross$-$sectional areas of the HSVs for the charging units

A$\_$du $=1$ # Effective cross$-$sectional areas of the HSVs for the discharging units

alpha $=1$ # Correction coefficient of the flow area of the valve port

sigma $=0\mathrm{.}5$ # Critical pressure ratio

d $=1$ # Duty cycles for the process

# PID Controller Parameters

Kp $=0\mathrm{.}1$

Ki $=0\mathrm{.}01$

Kd $=0\mathrm{.}01$

# Initial conditions

p$0=$ p$\_$atm

p$\_$t$0=$ p$\_$atm

Q$\_$m$1=0$

Q$\_$m$2=0$

Q$\_$mv $=0$

x$\_$vm $=0$

# Time points

t $=$ np$.$linspace$(0,10,1000)$

epsilon $=0\mathrm{.}000001$ # Small time offset to avoid division by zero at t$=0$

# Pressure Differential Equation

def model$($y$,$ t$)$:

p$,$ p$\_$t$,$ Q$\_$m$1,$ Q$\_$m$2,$ Q$\_$mv$,$ x$\_$vm $=$ y

f$1=($k $*$ Q$\_$m$1*$ R $*$ T$)/$ V

f$2=($k $*$ R $*$ T$)/$ V$\_$t

Q$\_$m $=$ Q$\_$m$1-$ Q$\_$m$2$

f$3=($p$\_$t $*$ S$\_$v$)/($R $*$ T$)*(1-$ p$\_$v $/$ p$\_$t$)$

t$\_$safe $=$ t $+$ epsilon # Use t$\_$safe instead of t to avoid division by zero

dpdt $=$ f$1-$ Kp $*($p $-$ p$0)-$ Ki $*($p $-$ p$0)*$ t$\_$safe $-$ Kd $*($p $-$ p$0)/$ t$\_$safe

dp$\_$tdt $=$ f$2*($Q$\_$m$2-$ Q$\_$mv$)$

dQ$\_$m$1$dt $=0\mathrm{.}0404*$ A$\_$cu $*$ p$\_$atm $/$ np$.$sqrt$($T$)*$ f$($p$,$ p$\_$atm$)*$ d$\_$cu$($t$)$

dQ$\_$m$2$dt $=0\mathrm{.}0404*$ A$\_$du $*$ p $/$ np$.$sqrt$($T$)*$ f$($p$,$ p$\_$t$)*$ d$\_$du$($t$)$

dQ$\_$mvdt $=$ f$3$

dx$\_$vmdt $=0\mathrm{.}0404*$ np$.$pi $*$ alpha $*$ D $*$ x$\_$vm $*$ p $/$ np$.$sqrt$($T$)*$ f$($p$,$ p$\_$t$)*$ d

return $[$dpdt$,$ dp$\_$tdt$,$ dQ$\_$m$1$dt$,$ dQ$\_$m$2$dt$,$ dQ$\_$mvdt$,$ dx$\_$vmdt$]$

# Function f$($p$1,$ p$2)$

def f$($p$1,$ p$2)$:

ratio $=$ p$2/($p$1*(1-$ sigma$))$

if ratio sigma:

return $1$

else:

adjusted$\_$ratio $=$ max$(1-$ ratio $**2,0)$ # Ensure the value under the square root is non$-$negative

return np$.$sqrt$($adjusted$\_$ratio$)$

# Duty cycles

def d$\_$cu$($t$)$:

return $0\mathrm{.}5$ # Duty cycles for the process of charging

def d$\_$du$($t$)$:

return $0\mathrm{.}5$ # Duty cycles for the process of discharging

# Solve ODE

y$0=[$p$0,$ p$\_$t$0,$ Q$\_$m$1,$ Q$\_$m$2,$ Q$\_$mv$,$ x$\_$vm$]$

y $=$ odeint$($model$,$ y$0,$ t$)$

# Plot results

plt$.$figure$($figsize$=(10,5))$

plt$.$plot$($t$,$ y$[$:$,0],\text{'}$r$\text{'},$ label$=\text{'}$p$\text{'})$

plt$.$plot$($t$,$ y$[$:$,1],\text{'}$b$\text{'},$ label$=\text{'}$p$\_$t$\text{'})$

plt$.$legend$($loc$=$'best'$)$

plt$.$xlabel$(\text{'}$Time $($s$)\text{'})$

plt$.$ylabel$(\text{'}$Pressure $($Pa$)\text{'})$

plt$.$title$(\text{'}$Pressure Control'$)$

plt$.$grid$()$

plt$.$show$()$