temperature.py

   1 import Melcor
   2 import time
   3 import data_logger
   4 import stripchart
   5
   6 log_dir = "/home/wking/rsrch/data/temperature"
   7
   8 # buzzwords: 'integrator windup' for integral term built up during a slow approach.
   9
  10
  11 class error (Exception) :
  12     "Errors with the temperature controller"
  13     pass
  14
  15 class errorMelcor (error) :
  16     pass
  17 class errorOutOfRange (error) :
  18     pass
  19
  20 def _check1(functionCall) :
  21     (err, val) = functionCall
  22     if err != 0 :
  23         raise errorMelcor
  24     return val
  25
  26 def _check0(functionCall) :
  27     err = functionCall
  28     if err != 0 :
  29         raise errorMelcor
  30
  31
  32 def melcor2double(value) :
  33     (err, doub) = Melcor.melcor2double(value)
  34     if err != 0 :
  35         raise errorMelcor, "Error converting melcor to double"
  36     return doub
  37 def double2melcor(doub) :
  38     (err, val) = Melcor.double2melcor(doub)
  39     if err != 0 :
  40         raise errorMelcor, "Error converting double to melcor"
  41     return val
  42
  43 def check_range(raw_output, min, max) :
  44     if raw_output < min :
  45         raise errorOutOfRange, '%g < %g' % (raw_output, min)
  46     if raw_output > max :
  47         raise errorOutOfRange, '%g > %g' % (raw_output, max)
  48
  49 class tempController :
  50     "Pretty wrappers for controlling a Melcor MTCA Temperature Controller"
  51     def __init__(self, controller=1, device='/dev/ttyS0', maxCurrent=0.2) :
  52         """
  53         (controller, device, maxCurrent) -> (tempController instance)
  54         controller : MTCA controller Id
  55         device     : serial port you're using to connect to the controller
  56         maxCurrent : initial maximum allowed current (in Amps)
  57         Set maxCurrent = None if you don't want to adjust from it's prev. value.
  58
  59         0.2 A is the default max current since it seems ok to use without fluid
  60         cooled heatsink.  If you are cooling the heatsink, use 1.0 A, which seems
  61         safely below the peltier's 1.2 A limit.
  62         """
  63         self.verbose = False
  64         self.setpoint = 20.0 # degrees C
  65         self.Tmin = 5.0 # setup some protective bounds for sanity checks
  66         self.Tmax = 50.0
  67         self.specMaxCur = 4.0 # Amps, the rated max current from controller specs
  68         self.T = Melcor.tempController(controller, device)
  69         self.Tstrip = stripchart.stripchart(pipename='Tstrip_pipe',
  70                                             title='Temp strip')
  71         self.Cstrip = stripchart.stripchart(pipename='Cstrip_pipe',
  72                                             title='Current strip')
  73         if maxCurrent != None : # if None, just leave maxCurrent at it's prev. val.
  74             self.setMaxCurrent(maxCurrent) # Amps
  75     def getTemp(self) :
  76         "Returns the current process temperature in degrees Celsius"
  77         val = self.read(Melcor.REG_HIGH_RESOLUTION)
  78         temp = val/100.0
  79         if self.Tstrip.status == 'open' :
  80             self.Tstrip.add_point(temp)
  81         return temp
  82     def getAmbientTemp(self) :
  83         "Returns room temperature in degrees Celsius"
  84         val = self.read(Melcor.REG_AMBIENT_TEMPERATURE)
  85         # convert (Fahrenheit*10) to Celsius
  86         return (val/10.0 - 32)/1.8
  87     def setSetpoint(self, setpoint) :
  88         "Set the temperature setpoint in degrees Celsius"
  89         val = double2melcor(setpoint)
  90         self.write(Melcor.REG_SET_POINT_1, val)
  91     def getSetpoint(self) :
  92         "Get the temperature setpoint in degrees Celsius"
  93         val = self.read(Melcor.REG_SET_POINT_1)
  94         return melcor2double(val)
  95     def setMaxCurrent(self, maxCur) :
  96         """
  97         Set the max current in Amps.
  98         (Note to Melcor enthusiasts: set's both the 'above' and 'below' limits)
  99         """
 100         maxPercent = maxCur / self.specMaxCur * 100
 101         val = double2melcor(maxPercent)
 102         self.write(Melcor.REG_HIGH_POWER_LIMIT_ABOVE, val)
 103         self.write(Melcor.REG_HIGH_POWER_LIMIT_BELOW, val)
 104         self.maxCurrent = maxCur
 105     def getMaxCurrent(self) :
 106         """
 107         () -> (currentLimitAbove (A), currentLimitBelow (A), currentLimitSetpoint (deg C))
 108         """
 109         per = self.read(Melcor.REG_HIGH_POWER_LIMIT_ABOVE)
 110         curLimAbove = melcor2double(per)/100.0 * self.specMaxCur
 111         per = self.read(Melcor.REG_HIGH_POWER_LIMIT_BELOW)
 112         curLimBelow = melcor2double(per)/100.0 * self.specMaxCur
 113         val = self.read(Melcor.REG_POWER_LIMIT_SETPOINT)
 114         curLimSet = melcor2double(val)
 115         return (curLimAbove, curLimBelow, curLimSet)
 116     def getPercentCurrent(self) :
 117         """
 118         Returns the percent of rated max current being output.
 119         See getCurrent()
 120         """
 121         val = self.read(Melcor.REG_PERCENT_OUTPUT)
 122         return val
 123     def getCurrent(self) :
 124         """
 125         The returned current is not the actual current,
 126         but the current that the temperature controller
 127         calculates it should generate.
 128         If the voltage required to generate that current
 129         exceeds the controllers max voltage (15V on mine),
 130         then the physical current will be less than the
 131         value returned here.
 132         """
 133         percentOutput = self.getPercentCurrent()
 134         cur = self.specMaxCur * percentOutput / 100.0
 135         if self.Cstrip.status == 'open' :
 136             self.Cstrip.add_point(cur)
 137         return cur
 138     def setCoolingGains(self, propband=0.1, integral=0, derivative=0) :
 139         """
 140         (propband, integral, derivative, dead_band) -> None
 141         propband   : propotional gain band in degrees C
 142         integral   : integral weight in minutes (0.00 to 99.99)
 143         derivative : derivative weight in minutes (? to ?)
 144         See 5.10 and the pages afterwards in the manual for Melcor's explaination.
 145         Formula (from Cornell BioPhys El Producto Beamline notes)
 146         P_cout = -1/T_prop * [ (T_samp - T_set)
 147                                + 1/t_int * int_-inf^t (T_samp(t')-T_set(t')) dt'
 148                                + t_deriv * dT_samp/dt
 149         Where P_cout is the percent of the rated max current that the controller
 150          would like to output if you weren't limiting it,
 151         T_prop is the propband input to this function,
 152         T_samp is the measured temperature of the sample in deg C,
 153         T_set is the setpoint in deg C,
 154         t_int is the integral input to this function,
 155         the integral with respect to t' is actually only from the time that
 156          T_samp has been with T_prop of T_set (not -inf), and
 157         t_deriv is the derivative input to this function.
 158
 159         Cooling is output 1
 160         """
 161         check_range(propband, 0, 99.9)
 162         check_range(integral, 0, 99.99)
 163         check_range(derivative, 0, 99.99)
 164
 165         val = double2melcor(propband)
 166         self.write(Melcor.REG_PROPBAND_1, val)
 167         val = int(integral * 100)
 168         self.write(Melcor.REG_INTEGRAL_1, val)
 169         val = int(derivative * 100)
 170         self.write(Melcor.REG_DERIVATIVE_1, val)
 171     def getCoolingGains(self) :
 172         "() -> (propband, integral, derivative)"
 173         val = self.read(Melcor.REG_PROPBAND_1)
 174         propband = melcor2double(val)
 175         val = self.read(Melcor.REG_INTEGRAL_1)
 176         integral = val/100.0
 177         val = self.read(Melcor.REG_DERIVATIVE_1)
 178         derivative = val/100.0
 179         return (propband, integral, derivative)
 180     def setHeatingGains(self, propband=0.1, integral=0, derivative=0) :
 181         """
 182         (propband, integral, derivative, dead_band) -> None
 183         propband   : propotional gain band in degrees C
 184         integral   : integral weight in minutes (0.00 to 99.99)
 185         derivative : derivative weight in minutes (? to ?)
 186         Don't use derivative, dead time.
 187         Cycle time?
 188         Histerysis?
 189         Burst?
 190         See 5.10 and the pages afterwards in the manual for Melcor's explaination.
 191         Formula (from Cornell BioPhys El Producto Beamline notes)
 192         P_cout = -1/T_prop * [ (T_samp - T_set)
 193                                + 1/t_int * int_-inf^t (T_samp(t')-T_set(t')) dt'
 194                                + t_deriv * dT_samp/dt
 195         Where P_cout is the percent of the rated max current that the controller
 196          would like to output if you weren't limiting it,
 197         T_prop is the propband input to this function,
 198         T_samp is the measured temperature of the sample in deg C,
 199         T_set is the setpoint in deg C,
 200         t_int is the integral input to this function,
 201         the integral with respect to t' is actually only from the time that
 202          T_samp has been with T_prop of T_set (not -inf), and
 203         t_deriv is the derivative input to this function.
 204
 205         Heating is output 2
 206         """
 207         check_range(propband, 0, 99.9)
 208         check_range(integral, 0, 99.99)
 209         check_range(derivative, 0, 99.99)
 210
 211         val = double2melcor(propband)
 212         self.write(Melcor.REG_PROPBAND_2, val)
 213         val = int(integral * 100)
 214         self.write(Melcor.REG_INTEGRAL_2, val)
 215         val = int(derivative * 100)
 216         self.write(Melcor.REG_DERIVATIVE_2, val)
 217     def getHeatingGains(self) :
 218         "() -> (propband, integral, derivative)"
 219         val = self.read(Melcor.REG_PROPBAND_2)
 220         propband = melcor2double(val)
 221         val = self.read(Melcor.REG_INTEGRAL_2)
 222         integral = val/100.0
 223         val = self.read(Melcor.REG_DERIVATIVE_2)
 224         derivative = val/100.0
 225         return (propband, integral, derivative)
 226     def getFeedbackTerms(self) :
 227         """
 228         Experimental
 229         """
 230         pid = self.read(Melcor.REG_PID_POWER_1)
 231         prop = self.read(Melcor.REG_PROP_TERM_1)
 232         ntgrl = self.read(Melcor.REG_INTEGRAL_TERM_1)
 233         deriv = self.read(Melcor.REG_DERIVATIVE_TERM_1)
 234         pout = self.getPercentCurrent()
 235         temp = self.getTemp()
 236         tset = self.getSetpoint()
 237         print "pid %g =? sum %g =? cur %g" % (pid, prop+ntgrl+deriv, pout)
 238         print "meas:     prop %d, integral %d, deriv %d" % (prop, ntgrl, deriv)
 239         print "my calcs: prop %d" % (temp-tset)
 240     def setTemp(self, setpoint, tolerance=0.3, time=10.0) :
 241         """
 242         Changes setpoint to SETPOINT and waits for stability
 243         """
 244         self.setSetpoint(setpoint)
 245         while self.isStable(setpoint, tolerance, time) != True :
 246             pass
 247     def setTemp_funkygain(self, setpoint, dead_time, heat_rate, cool_rate,
 248                           peltier_efficiency_fn, outside_equilib_rate,
 249                           tolerance=0.3, time=10.0) :
 250         """
 251         Highly experimental, see diffusion.py
 252         """
 253         mode = ""
 254         T = self.getTemp()
 255         # full steam ahead
 256         print "full steam ahead"
 257         self.setSetpoint(setpoint)
 258         self.setHeatingGains(0.1, 0, 0)
 259         self.setCoolingGains(0.1, 0, 0)
 260         if T < setpoint :
 261             mode = "Heating"
 262             self._heat_until_close(setpoint, dead_time, heat_rate)
 263         elif T > setpoint :
 264             mode = "Cooling"
 265             self._cool_until_close(setpoint, dead_time, cool_rate)
 266         # coast
 267         print "coast while temperature equilibrates"
 268         self.setHeatingGains(100, 0, 0)
 269         self.setCoolingGains(100, 0, 0)
 270         time.sleep(dead_time*2)
 271         cool_prop, heat_prop = self.calcPropBands()
 272         print "calculated prop bands: c %g, h %g deg C" % (cool_prop, heat_prop)
 273         print "reset integral gain, and bump to predicted props"
 274         # pop down to reset integral gain, could also jump setpoint...
 275         self.setHeatingGains(0.1, 0, 0)
 276         self.setHeatingGains(heat_prop, 0, 0)
 277         self.setCoolingGains(0.1, 0, 0)
 278         self.setCoolingGains(cool_prop, 0, 0)
 279         time.sleep(dead_time*4)
 280         # now add in some integral to reduce droop
 281         print "set integral gains to %g" % (dead_time*4)
 282         self.setHeatingGains(heat_prop, dead_time*4, 0)
 283         self.setCoolingGains(cool_prop, dead_time*4, 0)
 284         time.sleep(dead_time*8)
 285         print "wait to enter tolerance band"
 286         while (self.getTemp()-setpoint) :
 287             time.sleep(dead_time)
 288         print "should be stable now"
 289         if not self.isStable(setpoint, tolerance, time) :
 290             raise error, "Algorithm broken ;)"
 291     def _heat_until_close(self, setpoint, dead_time, heat_rate) :
 292         while self.getTemp() < setpoint - 0.5*rate*dead_time :
 293             time.sleep(dead_time/10.0)
 294     def calcPropBands(setpoint, peltier_efficiency_fn, outside_equilib_rate) :
 295         heat_loss = outside_equilib_rate * (setpoint - self.getAmbientTemp())
 296         required_current = heat_loss / peltier_efficiency_fn(setpoint)
 297         if required_current > self.maxCurrent :
 298             raise errorOutOfRange, "Can't source %g Amps", required_current
 299         fraction_current = required_current / self.maxCurrent
 300         droop = 0.5 # expected droop in deg C on only proporitional gain
 301         # droop / T_prop = fraction current
 302         T_prop = droop / fraction_current
 303         if setpoint > self.getAmbientTemp()+5 : # heating
 304             return (T_prop*10, T_prop)
 305         elif setpoint < self.getAmbientTemp()+5 : # cooling
 306             return (T_prop, T_prop*10)
 307         else : # right about room temperature
 308             return (T_prop, T_prop)
 309     def stripT(self, on) :
 310         if on :
 311             self.stripT = True
 312             pipename = 'Temp_pipe'
 313             self.stripTpipe = os.popen(pipename)
 314             os.system("stripchart -u %s" % pipename)
 315     def isStable(self, setpoint, tolerance=0.3, maxTime=10.0) :
 316         """
 317         Counts how long the temperature stays within
 318         TOLERANCE of SETPOINT.
 319         Returns when temp goes bad, or MAXTIME elapses.
 320         """
 321         stable = False
 322         startTime = time.time()
 323         stopTime = startTime
 324         while abs(self.getTemp() - setpoint) < tolerance :
 325                 stopTime = time.time()
 326                 if (stopTime-startTime) > maxTime :
 327                         print "Stable for long enough"
 328                         break
 329         if stopTime-startTime > maxTime :
 330                 return True
 331         else :
 332                 return False
 333     def setFilterTime(self, seconds) :
 334         """
 335         Positive values to affect only monitored values.
 336         Negative values affect both monitored and control values.
 337         """
 338         decSeconds = int(seconds*10)
 339         if decSeconds < 0 : # convert (unsigned int) -> (2's compliment signed)
 340                 decSeconds += 2**16
 341         self.write(Melcor.REG_INPUT_SOFTWARE_FILTER_1, decSeconds)
 342     def getFilterTime(self) :
 343         """
 344         Positive values to affect only monitored values.
 345         Negative values affect both monitored and control values.
 346         """
 347         val = self.read(Melcor.REG_INPUT_SOFTWARE_FILTER_1)
 348         if val >= 2**15 : # convert (2's complement signed) -> (unsigned int)
 349                 val -= 2**16
 350         return val/10.0
 351     def sanityCheck(self) :
 352         "Check that some key registers have the values we expect"
 353         self._sanityCheck(Melcor.REG_UNITS_TYPE,   2) # SI
 354         self._sanityCheck(Melcor.REG_C_OR_F,       1) # C
 355         self._sanityCheck(Melcor.REG_FAILURE_MODE, 2) # off
 356         self._sanityCheck(Melcor.REG_RAMPING_MODE, 0) # off
 357         self._sanityCheck(Melcor.REG_OUTPUT_1,     1) # cool
 358         self._sanityCheck(Melcor.REG_OUTPUT_2,     1) # heat
 359     def _sanityCheck(self, register, expected_value) :
 360         val = self.read(register)
 361         if val != expected_value :
 362                 print "Register %d, expected %d, was %d" % (register,
 363                                                     expected_value,
 364                                                     val)
 365                 raise error, "Controller settings error"
 366     def read(self, register) :
 367         """
 368         (register) -> (value)
 369         Returns the value of the specified memory register on the controller.
 370         Registers are defined in the Melcor module.
 371         See melcor_registers.h for a pointers on meanings and manual page nums.
 372         """
 373         (err, val) = self.T.read(register)
 374         if err != 0 :
 375                 raise errorMelcor
 376         return val
 377     def write(self, register, value) :
 378         """
 379         (register, value) -> None
 380         Sets the value of the specified memory register on the controller.
 381         Registers are defined in the Melcor module.
 382         See melcor_registers.h for a pointers on meanings and manual page nums.
 383         """
 384         err = self.T.write(register, value)
 385         if err != 0 :
 386                 raise errorMelcor
 387     def getDeadtimeData(self, num_oscillations=10, curHysteresis=0.8, log=True) :
 388         orig_heat_gains = self.getHeatingGains()
 389         orig_cool_gains = self.getCoolingGains()
 390         if self.verbose :
 391             print "Measuring dead time"
 392             print " go to bang-bang"
 393         self.setHeatingGains(0.1, 0, 0)
 394         self.setCoolingGains(0.1, 0, 0)
 395         def isHeating(cur) :
 396             if cur > curHysteresis :
 397                 return True
 398             elif cur < -curHysteresis :
 399                 return False
 400             else :
 401                 return None
 402         i=0
 403         timeArr = [0.0]
 404         temp = self.getTemp()
 405         cur = self.getCurrent()
 406         heat_first = isHeating(cur)
 407         start_time = time.time()
 408         tm = 0
 409         if verbose :
 410             print " Wait to exit hysteresis region"
 411         while heat_first == None and tm < 30:
 412                 temp = t.getTemp()
 413                 cur = t.getCurrent()
 414                 heat_first = isHeating(temp, cur)
 415                 tm = time.time()-start_time
 416         if tm > 30 :
 417                 raise error, "after 30 seconds, still inside hysteresis region"
 418         if self.verbose :
 419             print " Read oscillations"
 420         heating = heat_first
 421         start_time = time.time()
 422         tempArr = [temp]
 423         curArr = [cur]
 424         if verbose :
 425                 print "Temp %g\t(%g),\tCur %g,\tTime %d" % (temp, temp-Tset, cur, 0)
 426         while i < numOscillations*2 :
 427                 temp = t.getTemp()
 428                 tm = time.time()-start_time
 429                 cur = t.getCurrent()
 430                 tempArr.append(temp)
 431                 timeArr.append(tm)
 432                 curArr.append(cur)
 433                 check_signs(temp,cur)
 434                 if heating == True and isHeating(temp, cur) == False :
 435                         print "Transition to cooling (i=%d)" % i
 436                         heating = False
 437                         i += 1
 438                 elif heating == False and isHeating(temp, cur) == True :
 439                         print "Transition to heating (i=%d)" % i
 440                         heating = True
 441                         i += 1
 442         if verbose :
 443             print " Restoring gains"
 444         self.setHeatingGains(*orig_heat_gains)
 445         self.setCoolingGains(*orig_cool_gains)
 446         if log == True :
 447             log = 1# MARK
 448     def time_getTemp(self) :
 449         "Rough estimate timeing of getTemp(), takes me about 0.1s"
 450         start = time.time()
 451         for i in range(10) :
 452             self.getTemp()
 453         stop = time.time()
 454         return (stop-start)/10.0
 455
 456 def _test_tempController() :
 457     t = tempController(controller=1, maxCurrent=0.1)
 458
 459     print "Temp     = %g" % t.getTemp()
 460     print "Current  = %g" % t.getCurrent()
 461     print "Setpoint = %g" % t.getSetpoint()
 462
 463     print "Setting setpoint to 5.0 deg C"
 464     t.setSetpoint(5.0)
 465     sp = t.getSetpoint()
 466     print "Setpoint = %g" % sp
 467     if sp != 5.0 :
 468         raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0)
 469     time.sleep(10) # give the controller some time to overcome any integral gain
 470     c = t.getCurrent()
 471     print "Current  = %g" % c
 472     mca, mcb, mct = t.getMaxCurrent()
 473     if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb
 474         if c != mcb :
 475             raise Exception, "Current not at max %g, and we're shooting for a big temp" % mcb
 476     else :
 477         if c != mca :
 478             raise Exception, "Current not at max %g, and we're shooting for a big temp" % mca
 479
 480
 481     print "Setting setpoint to 50.0 deg C"
 482     t.setSetpoint(50.0)
 483     sp = t.getSetpoint()
 484     print "Setpoint = %g" % sp
 485     if sp != 5.0 :
 486         raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0)
 487     time.sleep(10)
 488     c = t.getCurrent()
 489     print "Current  = %g" % c
 490     print "Success"
 491     mca, mcb, mct = t.getMaxCurrent()
 492     if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb
 493         if -c != mcb :
 494             raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mcb)
 495     else :
 496         if -c != mca :
 497             raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mca)
 498
 499 def test() :
 500     _test_tempController()
 501
 502 if __name__ == "__main__" :
 503     test()