examples/pid_repsonse.py

   1 #!/usr/bin/env python
   2 # Copyright (C) 2011 W. Trevor King <wking@drexel.edu>
   3 #
   4 # This file is part of pypid.
   5 #
   6 # pypid is free software: you can redistribute it and/or
   7 # modify it under the terms of the GNU Lesser General Public
   8 # License as published by the Free Software Foundation, either
   9 # version 3 of the License, or (at your option) any later version.
  10 #
  11 # pypid is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU Lesser General Public License for more details.
  15 #
  16 # You should have received a copy of the GNU Lesser General Public
  17 # License along with pypid.  If not, see
  18 # <http://www.gnu.org/licenses/>.
  19
  20
  21 from argparse import ArgumentParser
  22 from sys import stdout
  23 from time import sleep
  24
  25 try:
  26     from matplotlib import pyplot
  27     from numpy import loadtxt
  28 except (ImportError,RuntimeError), e:
  29     pyplot = None
  30     loadtxt = None
  31     plot_import_error = e
  32
  33 from pypid.backend.test import TestBackend
  34 from pypid.rules import ziegler_nichols_step_response
  35
  36
  37 parser = ArgumentParser(description='Simulate a step response.')
  38 parser.add_argument(
  39     '-K', '--process-gain', metavar='GAIN', type=float, default=1,
  40     help='process gain (PV-units over MV-units)',)
  41 parser.add_argument(
  42     '-L', '--dead-time', metavar='TIME', type=float, default=1,
  43     help='system dead time (lag)')
  44 parser.add_argument(
  45     '-T', '--decay-time', metavar='TIME', type=float, default=1,
  46     help='exponential decay timescale')
  47 parser.add_argument(
  48     '-P', '--proportional', metavar='GAIN', type=float, default=None,
  49     help='process gain (output units over input units)',)
  50 parser.add_argument(
  51     '-I', '--integral', metavar='TIME', type=float, default=None,
  52     help='intergral gain timescale')
  53 parser.add_argument(
  54     '-D', '--derivative', metavar='TIME', type=float, default=None,
  55     help='derivative gain timescale')
  56 parser.add_argument(
  57     '-M', '--max-mv', metavar='MV', type=float, default=100.,
  58     help='maximum manipulated variable')
  59 parser.add_argument(
  60     '-A', '--tuning-algorithm', metavar='TUNER', default=None,
  61     choices=['ZN'], help='step tuning algorithm')
  62 parser.add_argument(
  63     '-m', '--mode', metavar='MODE', default='PID',
  64     choices=['P', 'PI', 'PID'], help='controller mode')
  65 parser.add_argument(
  66     '-t', '--time', metavar='TIME', type=float, default=10.,
  67     help='simulation time')
  68 parser.add_argument(
  69     '-o', '--output', default='-', help='output log file')
  70 parser.add_argument(
  71     '-p', '--plot', action='store_true', default=False,
  72     help='plot the repsonse')
  73
  74 args = parser.parse_args()
  75
  76 if args.plot and not (pyplot and loadtxt) :
  77     raise plot_import_error
  78
  79 if args.output == '-':
  80     log_stream = stdout
  81     if args.plot:
  82         raise ValueError('can only plot when outputing to a file')
  83 else:
  84     log_stream = open(args.output, 'w')
  85
  86 K = args.process_gain
  87 L = args.dead_time
  88 T = args.decay_time
  89
  90 p,i,d = (0, float('inf'), 0)
  91 if args.tuning_algorithm == 'ZN':
  92     p,i,d = ziegler_nichols_step_response(
  93         process_gain=K, dead_time=L, decay_time=T, mode=args.mode)
  94 else:
  95     if args.proportional:
  96         p = args.proportional
  97     if args.integral:
  98         i = args.integral
  99     if args.derivative:
 100         d = args.derivative
 101
 102 b = TestBackend(
 103     process_gain=K, dead_time=L, decay_time=T, max_mv=args.max_mv,
 104     log_stream=log_stream)
 105 try:
 106     b.set_up_gains(proportional=p, integral=i, derivative=d)
 107     b.set_down_gains(proportional=p, integral=i, derivative=d)
 108     b.set_setpoint(1.0)
 109     sleep(args.time)
 110 finally:
 111     b.cleanup();
 112     if args.output != '-':
 113         log_stream.close()
 114
 115 if args.plot:
 116     header = open(args.output, 'r').readline()
 117     label = header.strip('#\n').split('\t')
 118     data = loadtxt(args.output)
 119     times = data[:,0] - data[0,0]
 120     pyplot.hold(True)
 121     subplot = 1
 122     for i in range(1, len(label)):
 123         if i in [1, 4, 6]:
 124             if i:
 125                 pyplot.legend(loc='best')  # add legend to previous subplot
 126             pyplot.subplot(3, 1, subplot)
 127             subplot += 1
 128         pyplot.plot(times, data[:,i], '.', label=label[i])
 129     pyplot.legend(loc='best')
 130     pyplot.show()