temperature.py

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