7 # buzzwords: 'integrator windup' for integral term built up during a slow approach.
9 class error (Exception) :
10 "Errors with the temperature controller"
13 class errorMelcor (error) :
15 class errorOutOfRange (error) :
18 def _check1(functionCall) :
19 (err, val) = functionCall
24 def _check0(functionCall) :
30 def melcor2double(value) :
31 (err, doub) = Melcor.melcor2double(value)
33 raise errorMelcor, "Error converting melcor to double"
35 def double2melcor(doub) :
36 (err, val) = Melcor.double2melcor(doub)
38 raise errorMelcor, "Error converting double to melcor"
41 def check_range(raw_output, min, max) :
43 raise errorOutOfRange, '%g < %g' % (raw_output, min)
45 raise errorOutOfRange, '%g > %g' % (raw_output, max)
47 class tempController :
48 "Pretty wrappers for controlling a Melcor MTCA Temperature Controller"
49 def __init__(self, controller=1, device='/dev/ttyS0', maxCurrent=0.2) :
51 (controller, device, maxCurrent) -> (tempController instance)
52 controller : MTCA controller Id
53 device : serial port you're using to connect to the controller
54 maxCurrent : initial maximum allowed current (in Amps)
55 Set maxCurrent = None if you don't want to adjust from it's prev. value.
57 0.2 A is the default max current since it seems ok to use without fluid
58 cooled heatsink. If you are cooling the heatsink, use 1.0 A, which seems
59 safely below the peltier's 1.2 A limit.
62 self.setpoint = 20.0 # degrees C
63 self.Tmin = 5.0 # setup some protective bounds for sanity checks
65 self.specMaxCur = 4.0 # Amps, the rated max current from controller specs
66 self.T = Melcor.tempController(controller, device)
67 self.Tstrip = stripchart.stripchart(pipename='Tstrip_pipe',
69 self.Cstrip = stripchart.stripchart(pipename='Cstrip_pipe',
70 title='Current strip')
71 if maxCurrent != None : # if None, just leave maxCurrent at it's prev. val.
72 self.setMaxCurrent(maxCurrent) # Amps
74 "Returns the current process temperature in degrees Celsius"
75 val = self.read(Melcor.REG_HIGH_RESOLUTION)
77 if self.Tstrip.status == 'open' :
78 self.Tstrip.add_point(temp)
80 def getAmbientTemp(self) :
81 "Returns room temperature in degrees Celsius"
82 val = self.read(Melcor.REG_AMBIENT_TEMPERATURE)
83 # convert (Fahrenheit*10) to Celsius
84 return (val/10.0 - 32)/1.8
85 def setSetpoint(self, setpoint) :
86 "Set the temperature setpoint in degrees Celsius"
87 val = double2melcor(setpoint)
88 self.write(Melcor.REG_SET_POINT_1, val)
89 def getSetpoint(self) :
90 "Get the temperature setpoint in degrees Celsius"
91 val = self.read(Melcor.REG_SET_POINT_1)
92 return melcor2double(val)
93 def setMaxCurrent(self, maxCur) :
95 Set the max current in Amps.
96 (Note to Melcor enthusiasts: set's both the 'above' and 'below' limits)
98 maxPercent = maxCur / self.specMaxCur * 100
99 val = double2melcor(maxPercent)
100 self.write(Melcor.REG_HIGH_POWER_LIMIT_ABOVE, val)
101 self.write(Melcor.REG_HIGH_POWER_LIMIT_BELOW, val)
102 self.maxCurrent = maxCur
103 def getMaxCurrent(self) :
105 () -> (currentLimitAbove (A), currentLimitBelow (A), currentLimitSetpoint (deg C))
107 per = self.read(Melcor.REG_HIGH_POWER_LIMIT_ABOVE)
108 curLimAbove = melcor2double(per)/100.0 * self.specMaxCur
109 per = self.read(Melcor.REG_HIGH_POWER_LIMIT_BELOW)
110 curLimBelow = melcor2double(per)/100.0 * self.specMaxCur
111 val = self.read(Melcor.REG_POWER_LIMIT_SETPOINT)
112 curLimSet = melcor2double(val)
113 return (curLimAbove, curLimBelow, curLimSet)
114 def getPercentCurrent(self) :
116 Returns the percent of rated max current being output.
119 val = int(self.read(Melcor.REG_PERCENT_OUTPUT))
122 return float(val)/10.0
123 def getCurrent(self) :
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
133 percentOutput = self.getPercentCurrent()
134 cur = self.specMaxCur * percentOutput / 100.0
135 if self.Cstrip.status == 'open' :
136 self.Cstrip.add_point(cur)
138 def setCoolingGains(self, propband=0.1, integral=0, derivative=0) :
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.
161 check_range(propband, 0, 99.9)
162 check_range(integral, 0, 99.99)
163 check_range(derivative, 0, 9.99)
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)
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) :
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.
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.
207 check_range(propband, 0, 99.9)
208 check_range(integral, 0, 99.99)
209 check_range(derivative, 0, 9.99)
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)
223 val = self.read(Melcor.REG_DERIVATIVE_2)
224 derivative = val/100.0
225 return (propband, integral, derivative)
226 def getFeedbackTerms(self) :
230 pid = int(self.read(Melcor.REG_PID_POWER_1))
233 prop = int(self.read(Melcor.REG_PROP_TERM_1))
236 ntgrl = int(self.read(Melcor.REG_INTEGRAL_TERM_1))
240 deriv = int(self.read(Melcor.REG_DERIVATIVE_TERM_1))
243 return (pid, prop, ntgrl, deriv)
244 def checkFeedbackTerms(self) :
245 pid, prop, ntgrl, deriv = self.getFeedbackTerms()
246 pout = self.getPercentCurrent()
248 Tset = self.getSetpoint()
249 if T > Tset : # cooling
250 Tprop, tint, tderiv = self.getCoolingGains()
252 Tprop, tint, tderiv = self.getHeatingGains()
253 print "pid(read) %g =? sum(calc from terms) %g =? cur(read) %g" % (pid, prop+ntgrl+deriv, pout)
254 print "read: prop %d, integral %d, deriv %d" % (prop, ntgrl, deriv)
255 print "my calcs: prop %g" % ((Tset-T)/Tprop)
256 def setTemp(self, setpoint, tolerance=0.3, time=10.0) :
258 Changes setpoint to SETPOINT and waits for stability
260 self.setSetpoint(setpoint)
261 while self.isStable(setpoint, tolerance, time) != True :
263 def setTemp_funkygain(self, setpoint, dead_time, heat_rate, cool_rate,
264 peltier_efficiency_fn, outside_equilib_rate,
265 tolerance=0.3, time=10.0) :
267 Highly experimental, see diffusion.py
272 print "full steam ahead"
273 self.setSetpoint(setpoint)
274 self.setHeatingGains(0.1, 0, 0)
275 self.setCoolingGains(0.1, 0, 0)
278 self._heat_until_close(setpoint, dead_time, heat_rate)
281 self._cool_until_close(setpoint, dead_time, cool_rate)
283 print "coast while temperature equilibrates"
284 self.setHeatingGains(100, 0, 0)
285 self.setCoolingGains(100, 0, 0)
286 time.sleep(dead_time*2)
287 cool_prop, heat_prop = self.calcPropBands()
288 print "calculated prop bands: c %g, h %g deg C" % (cool_prop, heat_prop)
289 print "reset integral gain, and bump to predicted props"
290 # pop down to reset integral gain, could also jump setpoint...
291 self.setHeatingGains(0.1, 0, 0)
292 self.setHeatingGains(heat_prop, 0, 0)
293 self.setCoolingGains(0.1, 0, 0)
294 self.setCoolingGains(cool_prop, 0, 0)
295 time.sleep(dead_time*4)
296 # now add in some integral to reduce droop
297 print "set integral gains to %g" % (dead_time*4)
298 self.setHeatingGains(heat_prop, dead_time*4, 0)
299 self.setCoolingGains(cool_prop, dead_time*4, 0)
300 time.sleep(dead_time*8)
301 print "wait to enter tolerance band"
302 while (self.getTemp()-setpoint) :
303 time.sleep(dead_time)
304 print "should be stable now"
305 if not self.isStable(setpoint, tolerance, time) :
306 raise error, "Algorithm broken ;)"
307 def _heat_until_close(self, setpoint, dead_time, heat_rate) :
308 while self.getTemp() < setpoint - 0.5*rate*dead_time :
309 time.sleep(dead_time/10.0)
310 def calcPropBands_HACK(setpoint, peltier_efficiency_fn, outside_equilib_rate) :
311 heat_loss = outside_equilib_rate * (setpoint - self.getAmbientTemp())
312 required_current = heat_loss / peltier_efficiency_fn(setpoint)
313 if required_current > self.maxCurrent :
314 raise errorOutOfRange, "Can't source %g Amps", required_current
315 fraction_current = required_current / self.maxCurrent
316 droop = 0.5 # expected droop in deg C on only proporitional gain
317 # droop / T_prop = fraction current
318 T_prop = droop / fraction_current
319 if setpoint > self.getAmbientTemp()+5 : # heating
320 return (T_prop*10, T_prop)
321 elif setpoint < self.getAmbientTemp()+5 : # cooling
322 return (T_prop, T_prop*10)
323 else : # right about room temperature
324 return (T_prop, T_prop)
326 mcode = self.read(Melcor.REG_AUTO_MANUAL_OP_MODE)
332 raise error, "Unrecognized mode code %d" % mcode
333 def setMode(self, mode) :
336 elif mode == 'manual' :
339 raise error, "Unrecognized mode %s" % mode
340 self.write(Melcor.REG_AUTO_MANUAL_OP_MODE, mcode)
341 def getManualCurrent(self) :
342 val = int(self.read(Melcor.REG_MANUAL_SET_POINT))
343 if val > 2**15 : # convert to signed
345 pct = float(val)/10.0 # stored value is percent * 10
346 return self.specMaxCur * pct / 100.0
347 def setManualCurrent(self, amps) :
348 if amps > self.maxCurrent :
349 raise error, "Suggested current %g > max %g" % \
350 (amps, self.maxCurrent)
351 pct = amps / self.specMaxCur * 100.0
352 val = int(pct * 10.0)
353 if val < 0 : # convert to unsigned
355 self.write(Melcor.REG_MANUAL_SET_POINT, val)
356 def calcPropBands_ZN_get_step_response(self,
357 initial_current=None,
358 initial_wait_time=20.0,
360 response_wait_time=40.0,
363 Ziegler-Nichols tuning, using step response for input.
364 Process must be stable when calling this function.
366 if initial_current == None :
367 initial_current = self.getCurrent()
368 original_mode = self.getMode()
369 if original_mode == 'manual' :
370 original_current = self.getManualCurrent()
372 self.setMode('manual')
373 self.setManualCurrent(initial_current)
378 # get some stability data before stepping
379 while tm < start + initial_wait_time :
380 tarr.append(tm-start)
381 Tarr.append(self.getTemp())
384 self.setManualCurrent(initial_current + current_step)
386 while tm < start + initial_wait_time + response_wait_time :
387 Tnow = self.getTemp()
388 if Tnow != Tlast : # save us some trouble averaging later
389 tarr.append(tm-start)
390 Tarr.append(self.getTemp())
394 self.setMode(original_mode)
395 if original_mode == 'manual' :
396 self.setManualCurrent(original_current)
397 if plotVerbose == True :
398 from pylab import figure, plot, title, xlabel, ylabel
400 plot(tarr, Tarr, 'r.-')
403 title('Plant step response')
405 def calcPropBands_ZN_analyze_step_response(self, tarr, Tarr,
406 initial_wait_time=20,
409 Analyze the step response to determine dead time td,
410 and slope at point of inflection spoi.
413 while tarr[i] < initial_wait_time :
415 # find point of inflection (steepest slope)
419 return (Tarr[i+1] - Tarr[i])/(tarr[i+1]-tarr[i])
420 for i in range(istep, len(tarr)-1) :
421 if slope(i) > slope(ipoi) :
423 print "Max slope at t = %g, T = %g" % (tarr[ipoi], Tarr[ipoi])
426 # find the initial temperature
428 for i in range(istep) :
430 initialT /= float(istep)
431 print "Initial temperature %g" % initialT
432 deltaT = Tarr[ipoi] - initialT
433 rise_t_to_poi = spoi/deltaT
434 td = tarr[ipoi] - initial_wait_time - rise_t_to_poi
436 def calcPropBands_ZN_compute_terms(self, td, spoi) :
441 def calcPropBands_ZN(self, initial_current=None,
442 initial_wait_time=20.0,
444 response_wait_time=40.0,
447 tarr, Tarr = self.calcPropBands_ZN_get_step_response( \
448 initial_current, initial_wait_time,
449 current_step, response_wait_time,
451 td, spoi = self.calcPropBands_ZN_analyze_step_response( \
452 tarr, Tarr, initial_wait_time, textVerbose)
453 return self.calcPropBands_ZN_compute_terms(td, spoi)
454 def stripT(self, on) :
457 pipename = 'Temp_pipe'
458 self.stripTpipe = os.popen(pipename)
459 os.system("stripchart -u %s" % pipename)
460 def isStable(self, setpoint, tolerance=0.3, maxTime=10.0) :
462 Counts how long the temperature stays within
463 TOLERANCE of SETPOINT.
464 Returns when temp goes bad, or MAXTIME elapses.
467 startTime = time.time()
469 while abs(self.getTemp() - setpoint) < tolerance :
470 stopTime = time.time()
471 if (stopTime-startTime) > maxTime :
472 print "Stable for long enough"
474 if stopTime-startTime > maxTime :
478 def setFilterTime(self, seconds) :
480 Positive values to affect only monitored values.
481 Negative values affect both monitored and control values.
483 decSeconds = int(seconds*10)
484 if decSeconds < 0 : # convert (unsigned int) -> (2's compliment signed)
486 self.write(Melcor.REG_INPUT_SOFTWARE_FILTER_1, decSeconds)
487 def getFilterTime(self) :
489 Positive values to affect only monitored values.
490 Negative values affect both monitored and control values.
492 val = self.read(Melcor.REG_INPUT_SOFTWARE_FILTER_1)
493 if val >= 2**15 : # convert (2's complement signed) -> (unsigned int)
496 def sanityCheck(self) :
497 "Check that some key registers have the values we expect"
498 self._sanityCheck(Melcor.REG_UNITS_TYPE, 2) # SI
499 self._sanityCheck(Melcor.REG_C_OR_F, 1) # C
500 self._sanityCheck(Melcor.REG_FAILURE_MODE, 2) # off
501 self._sanityCheck(Melcor.REG_RAMPING_MODE, 0) # off
502 self._sanityCheck(Melcor.REG_OUTPUT_1, 1) # cool
503 self._sanityCheck(Melcor.REG_OUTPUT_2, 1) # heat
504 def _sanityCheck(self, register, expected_value) :
505 val = self.read(register)
506 if val != expected_value :
507 print "Register %d, expected %d, was %d" % (register,
510 raise error, "Controller settings error"
511 def read(self, register) :
513 (register) -> (value)
514 Returns the value of the specified memory register on the controller.
515 Registers are defined in the Melcor module.
516 See melcor_registers.h for a pointers on meanings and manual page nums.
518 (err, val) = self.T.read(register)
522 def write(self, register, value) :
524 (register, value) -> None
525 Sets the value of the specified memory register on the controller.
526 Registers are defined in the Melcor module.
527 See melcor_registers.h for a pointers on meanings and manual page nums.
529 err = self.T.write(register, value)
532 def getDeadtimeData(self, num_oscillations=10, curHysteresis=0.8) :
533 orig_heat_gains = self.getHeatingGains()
534 orig_cool_gains = self.getCoolingGains()
536 print "Measuring dead time"
537 print " go to bang-bang"
538 self.setHeatingGains(0.1, 0, 0)
539 self.setCoolingGains(0.1, 0, 0)
541 if cur > curHysteresis :
543 elif cur < -curHysteresis :
549 temp = self.getTemp()
550 cur = self.getCurrent()
551 heat_first = isHeating(cur)
552 start_time = time.time()
555 print " Wait to exit hysteresis region"
556 while heat_first == None and tm < 30:
559 heat_first = isHeating(temp, cur)
560 tm = time.time()-start_time
562 raise error, "after 30 seconds, still inside hysteresis region"
564 print " Read oscillations"
566 start_time = time.time()
570 print "Temp %g\t(%g),\tCur %g,\tTime %d" % (temp, temp-Tset, cur, 0)
571 while i < numOscillations*2 :
573 tm = time.time()-start_time
578 check_signs(temp,cur)
579 if heating == True and isHeating(temp, cur) == False :
580 print "Transition to cooling (i=%d)" % i
583 elif heating == False and isHeating(temp, cur) == True :
584 print "Transition to heating (i=%d)" % i
588 print " Restoring gains"
589 self.setHeatingGains(*orig_heat_gains)
590 self.setCoolingGains(*orig_cool_gains)
591 def time_getTemp(self) :
592 "Rough estimate timeing of getTemp(), takes me about 0.1s"
597 return (stop-start)/10.0
599 def _test_tempController() :
600 t = tempController(controller=1, maxCurrent=0.1)
602 print "Temp = %g" % t.getTemp()
603 print "Current = %g" % t.getCurrent()
604 print "Setpoint = %g" % t.getSetpoint()
606 print "Setting setpoint to 5.0 deg C"
609 print "Setpoint = %g" % sp
611 raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0)
612 time.sleep(10) # give the controller some time to overcome any integral gain
614 print "Current = %g" % c
615 mca, mcb, mct = t.getMaxCurrent()
616 if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb
618 raise Exception, "Current not at max %g, and we're shooting for a big temp" % mcb
621 raise Exception, "Current not at max %g, and we're shooting for a big temp" % mca
624 print "Setting setpoint to 50.0 deg C"
627 print "Setpoint = %g" % sp
629 raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0)
632 print "Current = %g" % c
634 mca, mcb, mct = t.getMaxCurrent()
635 if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb
637 raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mcb)
640 raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mca)
643 _test_tempController()
645 if __name__ == "__main__" :