Chamath/pid controller

Firmware/RaspberryPi/backend-pi/MQTTTransceiver.py

@@ -9,13 +9,11 @@
 
 
 class MQTTTransceiver:
-    PRESSURE_TOPIC = 'Ventilator/pressure'
-    FLOWRATE_TOPIC = 'Ventilator/flow_rate'
-    VOLUME_TOPIC = 'Ventilator/volume'
     CHART_DATA_TOPIC = 'Ventilator/chart_data'
     ACTUAL_TIDAL_VOLUME_TOPIC = 'Ventilator/vt'
     MINUTE_VOLUME_TOPIC = 'Ventilator/minute_volume'
     PIP_TOPIC = 'Ventilator/pip'
+    PRESSURE_DATA_TOPIC = "Ventilator/pressure_data"
 
     CALIB_FLOW_RATE_CONFIG_TOPIC = 'Config/calib_flow_rate'
     FIO2_CONFIG_TOPIC = 'Config/fio2'

Firmware/RaspberryPi/backend-pi/PID.py

@@ -0,0 +1,129 @@
+#!/usr/bin/python
+#
+# This file is part of IvPID.
+# Copyright (C) 2015 Ivmech Mechatronics Ltd. <[email protected]>
+#
+# IvPID is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# IvPID is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# title           :PID.py
+# description     :python pid controller
+# author          :Caner Durmusoglu
+# date            :20151218
+# version         :0.1
+# notes           :
+# python_version  :2.7
+# ==============================================================================
+
+"""Ivmech PID Controller is simple implementation of a Proportional-Integral-Derivative (PID) Controller in the Python Programming Language.
+More information about PID Controller: http://en.wikipedia.org/wiki/PID_controller
+"""
+import time
+
+class PID:
+    """PID Controller
+    """
+
+    def __init__(self, P=0.2, I=0.0, D=0.0, current_time=None):
+
+        self.Kp = P
+        self.Ki = I
+        self.Kd = D
+
+        self.sample_time = 0.00
+        self.current_time = current_time if current_time is not None else time.time()
+        self.last_time = self.current_time
+
+        self.clear()
+
+    def clear(self):
+        """Clears PID computations and coefficients"""
+        self.SetPoint = 0.0
+
+        self.PTerm = 0.0
+        self.ITerm = 0.0
+        self.DTerm = 0.0
+        self.last_error = 0.0
+
+        # Windup Guard
+        self.int_error = 0.0
+        self.windup_guard = 20.0
+
+        self.output = 0.0
+
+    def update(self, feedback_value, current_time=None):
+        """Calculates PID value for given reference feedback
+
+        .. math::
+            u(t) = K_p e(t) + K_i \int_{0}^{t} e(t)dt + K_d {de}/{dt}
+
+        .. figure:: images/pid_1.png
+           :align:   center
+
+           Test PID with Kp=1.2, Ki=1, Kd=0.001 (test_pid.py)
+
+        """
+        error = self.SetPoint - feedback_value
+
+        self.current_time = current_time if current_time is not None else time.time()
+        delta_time = self.current_time - self.last_time
+        delta_error = error - self.last_error
+
+        if (delta_time >= self.sample_time):
+            self.PTerm = self.Kp * error
+            self.ITerm += error * delta_time
+
+            if (self.ITerm < -self.windup_guard):
+                self.ITerm = -self.windup_guard
+            elif (self.ITerm > self.windup_guard):
+                self.ITerm = self.windup_guard
+
+            self.DTerm = 0.0
+            if delta_time > 0:
+                self.DTerm = delta_error / delta_time
+
+            # Remember last time and last error for next calculation
+            self.last_time = self.current_time
+            self.last_error = error
+
+            self.output = self.PTerm + (self.Ki * self.ITerm) + (self.Kd * self.DTerm)
+
+    def setKp(self, proportional_gain):
+        """Determines how aggressively the PID reacts to the current error with setting Proportional Gain"""
+        self.Kp = proportional_gain
+
+    def setKi(self, integral_gain):
+        """Determines how aggressively the PID reacts to the current error with setting Integral Gain"""
+        self.Ki = integral_gain
+
+    def setKd(self, derivative_gain):
+        """Determines how aggressively the PID reacts to the current error with setting Derivative Gain"""
+        self.Kd = derivative_gain
+
+    def setWindup(self, windup):
+        """Integral windup, also known as integrator windup or reset windup,
+        refers to the situation in a PID feedback controller where
+        a large change in setpoint occurs (say a positive change)
+        and the integral terms accumulates a significant error
+        during the rise (windup), thus overshooting and continuing
+        to increase as this accumulated error is unwound
+        (offset by errors in the other direction).
+        The specific problem is the excess overshooting.
+        """
+        self.windup_guard = windup
+
+    def setSampleTime(self, sample_time):
+        """PID that should be updated at a regular interval.
+        Based on a pre-determined sampe time, the PID decides if it should compute or return immediately.
+        """
+        self.sample_time = sample_time

Firmware/RaspberryPi/backend-pi/RnD/PIDControllerTuner.py

@@ -0,0 +1,95 @@
+import sys
+sys.path.append('..')
+
+from PID import PID
+from SensorReaderService import SensorReaderService
+from Variables import Variables
+import time
+from os import path
+import RPi.GPIO as GPIO
+
+# Constants
+PWM_FREQ = 2  # frequency for PWM
+SO_PIN = 6   # PIN (PWM) for O2 intake solenoid
+SI_PIN = 12  # PIN (PWM) for inspiratory solenoid
+SE_PIN = 14  # PIN (PWM) for expiratory solenoid
+PWM_O = None
+PWM_I = None
+PWM_E = None
+pid = None
+
+
+def load_pid_config():
+    global pid
+    with open('./pid.conf', 'r') as f:
+        config = f.readline().split(',')
+        pid.SetPoint = float(Variables.ps)
+        pid.setKp(float(config[0]))
+        pid.setKi(float(config[1]))
+        pid.setKd(float(config[2]))
+
+
+def create_pid_config():
+    if not path.isfile('./pid.conf'):
+        with open('./pid.conf', 'w') as f:
+            f.write('%s,%s,%s' % (Variables.Kp, Variables.Ki, Variables.Kd))
+
+
+def convert_pressure(p_hpa):
+    """ returns inspiratory pressure relative to atm in cmH2O"""
+    return (p_hpa * 1.0197442) - 1033.23
+
+
+def init_parameters():
+    global PWM_O, PWM_I, PWM_E, pid
+
+    # Initialize PWM pins
+    GPIO.setmode(GPIO.BCM)
+    GPIO.setwarnings(False)
+    GPIO.setup(SO_PIN, GPIO.OUT)
+    GPIO.setup(SI_PIN, GPIO.OUT)
+    GPIO.setup(SE_PIN, GPIO.OUT)
+
+    PWM_O = GPIO.PWM(SO_PIN, PWM_FREQ)
+    PWM_I = GPIO.PWM(SI_PIN, PWM_FREQ)
+    PWM_E = GPIO.PWM(SE_PIN, PWM_FREQ)
+
+    # Start the sensor reading service
+    sensing_service = SensorReaderService()
+
+    # Initialize PID controller
+    create_pid_config()
+    pid = PID(Variables.Kp, Variables.Ki, Variables.Kd)
+    pid.SetPoint = Variables.ps
+    pid.setSampleTime(Variables.pid_sampling_period)
+
+    # Open all values
+    PWM_O.ChangeDutyCycle(100)
+    PWM_I.ChangeDutyCycle(100)
+    PWM_E.ChangeDutyCycle(0)    # Normally open, hence duty_ratio=0
+
+###################################################################
+
+
+init_parameters()
+
+while True:
+    # load the latest PID related config. [Kp, Ki, Kd]
+    load_pid_config()
+
+    # read pressure data
+    pressure = Variables.p3
+
+    if pressure is None:
+        continue
+
+    pid.update(convert_pressure(pressure))
+    target_duty_ratio = pid.output
+    target_duty_ratio = max(min(int(target_duty_ratio), 100), 0)
+
+    print("Target: %.1f | Current: %.1f | Duty Ratio: %d" % (Variables.ps, convert_pressure(pressure), target_duty_ratio))
+
+    # Set PWM to target duty
+    PWM_I.ChangeDutyCycle(target_duty_ratio)
+
+    time.sleep(Variables.pid_sampling_period - 0.005)

Firmware/RaspberryPi/backend-pi/RnD/pid.conf

@@ -0,0 +1 @@
+10,1,1

Firmware/RaspberryPi/backend-pi/SensorReader.py

@@ -2,7 +2,6 @@
 import time
 import bme680
 from Variables import Variables
-from datetime import datetime
 
 
 class SensorReader:
@@ -64,7 +63,9 @@ def read_temp(self):
         self.fTemp = self.cTemp * 1.8 + 32
 
     def read_pressure(self):
-        if (self.bus_number == Variables.BUS_3):
+
+        # For demo p3 is read from a bme680 sensor
+        if Variables.demo and self.bus_number == Variables.BUS_3:
             return self.read_bme680()
 
         #Reading Data from i2c bus 3

Firmware/RaspberryPi/backend-pi/Variables.py

@@ -5,18 +5,32 @@
 
 
 class Variables:
+    # Pressure sensor BUS ids
     BUS_1 = 1 # first inspiratory pressure sensor
     BUS_2 = 3 # second inspiratory pressure sensor
     BUS_3 = 4 # first expiratory pressure sensor
     BUS_4 = 5 # second expiratory pressure sensor
-    calib_flow_rate = 10    # flow rate for calibration
+
+    # Pressure sensor values
+    p1 = 0      # pressure sensor 1 on inspiratory phase
+    p2 = 0      # pressure sensor 2 on inspiratory phase
+    p3 = 0      # pressure sensor 3 on expiratory phase
+    p4 = 0      # pressure sensor 4 on expiratory phase
+
+    # From GUI
+    calib_flow_rate = -1    # flow rate used for calibration. Set by GUI
     fio2 = 0    # fio2 value
     ie = 1      # I:E ratio
     rr = 10     # respiratory rate (RR)
     vt = 500    # tidal volume (mL)
     peep = 10   # PEEP
     ps = 20     # pressure support (Ps)
-    p1 = 0      # pressure sensor 1 on inspiratory phase
-    p2 = 0      # pressure sensor 2 on inspiratory phase
-    p3 = 0      # pressure sensor 3 on expiratory phase
-    p4 = 0      # pressure sensor 4 on expiratory phase
+
+    # PID controller parameters
+    Kp = 10
+    Ki = 1
+    Kd = 1
+    pid_sampling_period = 0.2
+
+    # Flag to indicate demo setup, where P3 is from a BME680 sensor
+    demo = False

Firmware/RaspberryPi/frontend-pi/public/index.html

@@ -130,13 +130,13 @@
         <div class="modal-dialog" role="document">
           <div class="modal-content">
             <div class="modal-header">
-              <h5 class="modal-title" id="calibarateModalLabel">Flow control calibaration</h5>
+              <h5 class="modal-title" id="calibarateModalLabel">Flow control calibration</h5>
               <button type="button" class="close" data-dismiss="modal" aria-label="Close">
                 <span aria-hidden="true">&times;</span>
               </button>
             </div>
             <div class="modal-body">
-              <p>Please enter flow control value to calibarate the air flow.</p>
+              <p>Please enter flow rate in L/min.</p>
               <input type="tel" class="keyboard form-control" id="calibarate_flow_input" type="number" onclick="show_easy_numpad(this);">
             </div>
             <div class="modal-footer">
@@ -257,7 +257,7 @@ <h5 class="modal-title" id="calibarateModalLabel">Flow control calibaration</h5>
                   <div class="card bg-dark text-center">
                     <div class="card-body">
                       <div>IE Ratio</div>
-                      <div id="ie_ratio" style="font-size:60px;font-style:initial;">1:2</div>
+                      <div id="ie_ratio" style="font-size:60px;font-style:initial;">1:1</div>
                     </div>
                   </div>
                 </td>
@@ -403,7 +403,9 @@ <h5 class="modal-title" id="calibarateModalLabel">Flow control calibaration</h5>
         },
         xAxis: {
           type: 'datetime',
+          tickPixelInterval: 150,
           dateTimeLabelFormats: {
+            millisecond: '%S.%L',
             second: '%S'
           },
         },

Firmware/RaspberryPi/frontend-pi/routes/tableValues.js

@@ -5,19 +5,7 @@ var router = express.Router();
 
 router.get('/values', function(req, res) {
 
-  //get Minute Volume by getting sum of last minute delivered volume
-  volumeData = database.get_volume();
-  pressure = database.get_pressure();
-  min_volume = 0;
-  for(var i = 0; i < volumeData.length; i++) {
-          min_volume = min_volume + volumeData[i][1];
-  }
-  max_pressure = 0;
-  for(var i = 0; i < pressure.length; i++) {
-          if (max_pressure < pressure[i][1]){
-            max_pressure = pressure[i][1];
-          }
-  }
+  //return the latest values
   res.json({
     "vt": database.get_vt(),
     "minute_volume": database.get_minute_volume(),

Firmware/RaspberryPi/frontend-pi/utils/mqtt_messenger.js

@@ -16,8 +16,10 @@ var self = module.exports = {
   mqtt_receiver: function() {
     console.log("Receiver called");
     client.on('connect', () => {
-      client.subscribe(ACTUAL_TIDAL_VOLUME_TOPIC)
       client.subscribe(CHART_DATA_TOPIC)
+      client.subscribe(ACTUAL_TIDAL_VOLUME_TOPIC)
+      client.subscribe(MINUTE_VOLUME_TOPIC)
+      client.subscribe(PIP_TOPIC)
     });
 
     client.on('message', (topic, message) => {
@@ -50,18 +52,18 @@ var self = module.exports = {
   },
 
   mqtt_vt: function(message) {
-    database.set_vt((parseFloat(message)).toFixed(2))
-    console.log("MQTT VT: " + (parseFloat(message)).toFixed(2));
+    database.set_vt(parseInt(message))
+    console.log("MQTT VT: " + parseInt(message));
   },
 
   mqtt_minute_volume: function(message) {
-    database.set_minute_volume((parseFloat(message)).toFixed(2))
-    console.log("MQTT Minute Volume: " + (parseFloat(message)).toFixed(2));
+    database.set_minute_volume(parseInt(message))
+    console.log("MQTT Minute Volume: " + parseInt(message));
   },
 
   mqtt_pip: function(message) {
-    database.set_pip((parseFloat(message)).toFixed(2))
-    console.log("MQTT Pip: " + (parseFloat(message)).toFixed(2));
+    database.set_pip((parseFloat(message)).toFixed(1))
+    console.log("MQTT Pip: " + (parseInt(message)));
   }
 
 };