From 0a6dd765f976634d92507a7916ab3aecf591434d Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Fri, 24 Oct 2008 07:31:43 -0400 Subject: [PATCH] Began versioning. --- temperature.py | 477 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 477 insertions(+) create mode 100644 temperature.py diff --git a/temperature.py b/temperature.py new file mode 100644 index 0000000..75c700c --- /dev/null +++ b/temperature.py @@ -0,0 +1,477 @@ +import Melcor +import time +import data_logger + +log_dir = "/home/wking/rsrch/data/temperature" + +class error (Exception) : + "Errors with the temperature controller" + pass + +class errorMmelcor (error) : + pass +class errorOutOfRange (error) : + pass + +def _check1(functionCall) : + (err, val) = functionCall + if err != 0 : + raise errorMelcor + return val + +def _check0(functionCall) : + err = functionCall + if err != 0 : + raise errorMelcor + + +def melcor2double(value) : + (err, doub) = Melcor.melcor2double(value) + if err != 0 : + raise errorMelcor, "Error converting melcor to double" + return doub +def double2melcor(doub) : + (err, val) = Melcor.double2melcor(doub) + if err != 0 : + raise errorMelcor, "Error converting double to melcor" + return val + +def check_range(raw_output, min, max) : + if raw_output < min : + raise errorOutOfRange, '%g < %g' % (raw_output, min) + if raw_output > max : + raise errorOutOfRange, '%g > %g' % (raw_output, max) + +class tempController : + "Pretty wrappers for controlling a Melcor MTCA Temperature Controller" + def __init__(self, controller=1, device='/dev/ttyS0', maxCurrent=0.2) : + """ + (controller, device, maxCurrent) -> (tempController instance) + controller : MTCA controller Id + device : serial port you're using to connect to the controller + maxCurrent : initial maximum allowed current (in Amps) + Set maxCurrent = None if you don't want to adjust from it's prev. value. + + 0.2 A is the default max current since it seems ok to use without fluid + cooled heatsink. If you are cooling the heatsink, use 1.0 A, which seems + safely below the peltier's 1.2 A limit. + """ + self.verbose = False + self.setpoint = 20.0 # degrees C + self.Tmin = 5.0 # setup some protective bounds for sanity checks + self.Tmax = 50.0 + self.specMaxCur = 4.0 # Amps, the rated max current from controller specs + self.T = Melcor.tempController(controller, device) + if maxCurrent != None : # if None, just leave maxCurrent at it's prev. val. + self.setMaxCurrent(maxCurrent) # Amps + def getTemp(self) : + "Returns the current process temperature in degrees Celsius" + val = self.read(Melcor.REG_HIGH_RESOLUTION) + temp = val/100.0 + return temp + def getAmbientTemp(self) : + "Returns room temperature in degrees Celsius" + val = self.read(Melcor.REG_AMBIENT_TEMPERATURE) + # convert (Fahrenheit*10) to Celsius + return (val/10.0 - 32)/1.8 + def setSetpoint(self, setpoint) : + "Set the temperature setpoint in degrees Celsius" + val = double2melcor(setpoint) + self.write(Melcor.REG_SET_POINT_1, val) + def getSetpoint(self) : + "Get the temperature setpoint in degrees Celsius" + val = self.read(Melcor.REG_SET_POINT_1) + return melcor2double(val) + def setMaxCurrent(self, maxCur) : + """ + Set the max current in Amps. + (Note to Melcor enthusiasts: set's both the 'above' and 'below' limits) + """ + maxPercent = maxCur / self.specMaxCur * 100 + val = double2melcor(maxPercent) + self.write(Melcor.REG_HIGH_POWER_LIMIT_ABOVE, val) + self.write(Melcor.REG_HIGH_POWER_LIMIT_BELOW, val) + self.maxCurrent = maxCur + def getMaxCurrent(self) : + """ + () -> (currentLimitAbove (A), currentLimitBelow (A), currentLimitSetpoint (deg C)) + """ + per = self.read(Melcor.REG_HIGH_POWER_LIMIT_ABOVE) + curLimAbove = melcor2double(per)/100.0 * self.specMaxCur + per = self.read(Melcor.REG_HIGH_POWER_LIMIT_BELOW) + curLimBelow = melcor2double(per)/100.0 * self.specMaxCur + val = self.read(Melcor.REG_POWER_LIMIT_SETPOINT) + curLimSet = melcor2double(val) + return (curLimAbove, curLimBelow, curLimSet) + def getPercentCurrent(self) : + """ + Returns the percent of rated max current being output. + See getCurrent() + """ + val = self.read(Melcor.REG_PERCENT_OUTPUT) + return val + def getCurrent(self) : + """ + The returned current is not the actual current, + but the current that the temperature controller + calculates it should generate. + If the voltage required to generate that current + exceeds the controllers max voltage (15V on mine), + then the physical current will be less than the + value returned here. + """ + percentOutput = self.getPercentCurrent() + return self.specMaxCur * percentOutput / 100 + def setCoolingGains(self, propband=0.1, integral=0, derivative=0) : + """ + (propband, integral, derivative, dead_band) -> None + propband : propotional gain band in degrees C + integral : integral weight in minutes (0.00 to 99.99) + derivative : derivative weight in minutes (? to ?) + See 5.10 and the pages afterwards in the manual for Melcor's explaination. + Formula (from Cornell BioPhys El Producto Beamline notes) + P_cout = -1/T_prop * [ (T_samp - T_set) + + 1/t_int * int_-inf^t (T_samp(t')-T_set(t')) dt' + + t_deriv * dT_samp/dt + Where P_cout is the percent of the rated max current that the controller + would like to output if you weren't limiting it, + T_prop is the propband input to this function, + T_samp is the measured temperature of the sample in deg C, + T_set is the setpoint in deg C, + t_int is the integral input to this function, + the integral with respect to t' is actually only from the time that + T_samp has been with T_prop of T_set (not -inf), and + t_deriv is the derivative input to this function. + + Cooling is output 1 + """ + check_range(propband, 0, 99.9) + check_range(integral, 0, 99.99) + check_range(derivative, 0, 99.99) + + val = double2melcor(propband) + self.write(Melcor.REG_PROPBAND_1, val) + val = int(integral * 100) + self.write(Melcor.REG_INTEGRAL_1, val) + val = int(derivative * 100) + self.write(Melcor.REG_DERIVATIVE_1, val) + def getCoolingGains(self) : + "() -> (propband, integral, derivative)" + val = self.read(Melcor.REG_PROPBAND_1) + propband = melcor2double(val) + val = self.read(Melcor.REG_INTEGRAL_1) + integral = val/100.0 + val = self.read(Melcor.REG_DERIVATIVE_1) + derivative = val/100.0 + return (propband, integral, derivative) + def setHeatingGains(self, propband=0.1, integral=0, derivative=0) : + """ + (propband, integral, derivative, dead_band) -> None + propband : propotional gain band in degrees C + integral : integral weight in minutes (0.00 to 99.99) + derivative : derivative weight in minutes (? to ?) + Don't use derivative, dead time. + Cycle time? + Histerysis? + Burst? + See 5.10 and the pages afterwards in the manual for Melcor's explaination. + Formula (from Cornell BioPhys El Producto Beamline notes) + P_cout = -1/T_prop * [ (T_samp - T_set) + + 1/t_int * int_-inf^t (T_samp(t')-T_set(t')) dt' + + t_deriv * dT_samp/dt + Where P_cout is the percent of the rated max current that the controller + would like to output if you weren't limiting it, + T_prop is the propband input to this function, + T_samp is the measured temperature of the sample in deg C, + T_set is the setpoint in deg C, + t_int is the integral input to this function, + the integral with respect to t' is actually only from the time that + T_samp has been with T_prop of T_set (not -inf), and + t_deriv is the derivative input to this function. + + Heating is output 2 + """ + check_range(propband, 0, 99.9) + check_range(integral, 0, 99.99) + check_range(derivative, 0, 99.99) + + val = double2melcor(propband) + self.write(Melcor.REG_PROPBAND_2, val) + val = int(integral * 100) + self.write(Melcor.REG_INTEGRAL_2, val) + val = int(derivative * 100) + self.write(Melcor.REG_DERIVATIVE_2, val) + def getHeatingGains(self) : + "() -> (propband, integral, derivative)" + val = self.read(Melcor.REG_PROPBAND_2) + propband = melcor2double(val) + val = self.read(Melcor.REG_INTEGRAL_2) + integral = val/100.0 + val = self.read(Melcor.REG_DERIVATIVE_2) + derivative = val/100.0 + return (propband, integral, derivative) + def getFeedbackTerms(self) : + """ + Experimental + """ + pid = self.read(Melcor.REG_PID_POWER_1) + prop = self.read(Melcor.REG_PROP_TERM_1) + ntgrl = self.read(Melcor.REG_INTEGRAL_TERM_1) + deriv = self.read(Melcor.REG_DERIVATIVE_TERM_1) + pout = self.getPercentCurrent() + temp = self.getTemp() + tset = self.getSetpoint() + print "pid %g =? sum %g =? cur %g" % (pid, prop+ntgrl+deriv, pout) + print "meas: prop %d, integral %d, deriv %d" % (prop, ntgrl, deriv) + print "my calcs: prop %d" % (temp-tset) + def setTemp(self, setpoint, tolerance=0.3, time=10.0) : + """ + Changes setpoint to SETPOINT and waits for stability + """ + self.setSetpoint(setpoint) + while self.isStable(setpoint, tolerance, time) != True : + pass + def setTemp_funkygain(self, setpoint, dead_time, heat_rate, cool_rate, + peltier_efficiency_fn, outside_equilib_rate, + tolerance=0.3, time=10.0) : + """ + Highly experimental, see diffusion.py + """ + mode = "" + T = self.getTemp() + # full steam ahead + print "full steam ahead" + self.setSetpoint(setpoint) + self.setHeatingGains(0.1, 0, 0) + self.setCoolingGains(0.1, 0, 0) + if T < setpoint : + mode = "Heating" + self._heat_until_close(setpoint, dead_time, heat_rate) + elif T > setpoint : + mode = "Cooling" + self._cool_until_close(setpoint, dead_time, cool_rate) + # coast + print "coast while temperature equilibrates" + self.setHeatingGains(100, 0, 0) + self.setCoolingGains(100, 0, 0) + time.sleep(dead_time*2) + cool_prop, heat_prop = self.calcPropBands() + print "calculated prop bands: c %g, h %g deg C" % (cool_prop, heat_prop) + print "reset integral gain, and bump to predicted props" + # pop down to reset integral gain, could also jump setpoint... + self.setHeatingGains(0.1, 0, 0) + self.setHeatingGains(heat_prop, 0, 0) + self.setCoolingGains(0.1, 0, 0) + self.setCoolingGains(cool_prop, 0, 0) + time.sleep(dead_time*4) + # now add in some integral to reduce droop + print "set integral gains to %g" % (dead_time*4) + self.setHeatingGains(heat_prop, dead_time*4, 0) + self.setCoolingGains(cool_prop, dead_time*4, 0) + time.sleep(dead_time*8) + print "wait to enter tolerance band" + while (self.getTemp()-setpoint) : + time.sleep(dead_time) + print "should be stable now" + if not self.isStable(setpoint, tolerance, time) : + raise error, "Algorithm broken ;)" + def _heat_until_close(self, setpoint, dead_time, heat_rate) : + while self.getTemp() < setpoint - 0.5*rate*dead_time : + time.sleep(dead_time/10.0) + def calcPropBands(setpoint, peltier_efficiency_fn, outside_equilib_rate) : + heat_loss = outside_equilib_rate * (setpoint - self.getAmbientTemp()) + required_current = heat_loss / peltier_efficiency_fn(setpoint) + if required_current > self.maxCurrent : + raise errorOutOfRange, "Can't source %g Amps", required_current + fraction_current = required_current / self.maxCurrent + droop = 0.5 # expected droop in deg C on only proporitional gain + # droop / T_prop = fraction current + T_prop = droop / fraction_current + if setpoint > self.getAmbientTemp()+5 : # heating + return (T_prop*10, T_prop) + elif setpoint < self.getAmbientTemp()+5 : # cooling + return (T_prop, T_prop*10) + else : # right about room temperature + return (T_prop, T_prop) + def isStable(self, setpoint, tolerance=0.3, maxTime=10.0) : + """ + Counts how long the temperature stays within + TOLERANCE of SETPOINT. + Returns when temp goes bad, or MAXTIME elapses. + """ + stable = False + startTime = time.time() + stopTime = startTime + while abs(self.getTemp() - setpoint) < tolerance : + stopTime = time.time() + if (stopTime-startTime) > maxTime : + print "Stable for long enough" + break + if stopTime-startTime > maxTime : + return True + else : + return False + def setFilterTime(self, seconds) : + """ + Positive values to affect only monitored values. + Negative values affect both monitored and control values. + """ + decSeconds = int(seconds*10) + if decSeconds < 0 : # convert (unsigned int) -> (2's compliment signed) + decSeconds += 2**16 + self.write(Melcor.REG_INPUT_SOFTWARE_FILTER_1, decSeconds) + def getFilterTime(self) : + """ + Positive values to affect only monitored values. + Negative values affect both monitored and control values. + """ + val = self.read(Melcor.REG_INPUT_SOFTWARE_FILTER_1) + if val >= 2**15 : # convert (2's complement signed) -> (unsigned int) + val -= 2**16 + return val/10.0 + def sanityCheck(self) : + "Check that some key registers have the values we expect" + _sanityCheck(Melcor.REG_UNITS_TYPE, 2) # SI + _sanityCheck(Melcor.REG_C_OR_F, 1) # C + _sanityCheck(Melcor.REG_FAILURE_MODE, 2) # off + _sanityCheck(Melcor.REG_RAMPING_MODE, 0) # off + _sanityCheck(Melcor.REG_OUTPUT_1, 1) # cool + _sanityCheck(Melcor.REG_OUTPUT_2, 1) # heat + def _sanityCheck(self, register, expected_value) : + val = self.read(register) + if val != expected_value : + print "Register %d, expected %d, was %d" % (register, + expected_value, + val) + raise error, "Controller settings error" + def read(self, register) : + """ + (register) -> (value) + Returns the value of the specified memory register on the controller. + Registers are defined in the Melcor module. + See melcor_registers.h for a pointers on meanings and manual page nums. + """ + (err, val) = self.T.read(register) + if err != 0 : + raise errorMelcor + return val + def write(self, register, value) : + """ + (register, value) -> None + Sets the value of the specified memory register on the controller. + Registers are defined in the Melcor module. + See melcor_registers.h for a pointers on meanings and manual page nums. + """ + err = self.T.write(register, value) + if err != 0 : + raise errorMelcor + def getDeadtimeData(self, num_oscillations=10, curHysteresis=0.8, log=True) : + orig_heat_gains = self.getHeatingGains() + orig_cool_gains = self.getCoolingGains() + if self.verbose : + print "Measuring dead time" + print " go to bang-bang" + self.setHeatingGains(0.1, 0, 0) + self.setCoolingGains(0.1, 0, 0) + def isHeating(cur) : + if cur > curHysteresis : + return True + elif cur < -curHysteresis : + return False + else : + return None + i=0 + timeArr = [0.0] + temp = self.getTemp() + cur = self.getCurrent() + heat_first = isHeating(cur) + start_time = time.time() + tm = 0 + if verbose : + print " Wait to exit hysteresis region" + while heat_first == None and tm < 30: + temp = t.getTemp() + cur = t.getCurrent() + heat_first = isHeating(temp, cur) + tm = time.time()-start_time + if tm > 30 : + raise error, "after 30 seconds, still inside hysteresis region" + if self.verbose : + print " Read oscillations" + heating = heat_first + start_time = time.time() + tempArr = [temp] + curArr = [cur] + if verbose : + print "Temp %g\t(%g),\tCur %g,\tTime %d" % (temp, temp-Tset, cur, 0) + while i < numOscillations*2 : + temp = t.getTemp() + tm = time.time()-start_time + cur = t.getCurrent() + tempArr.append(temp) + timeArr.append(tm) + curArr.append(cur) + check_signs(temp,cur) + if heating == True and isHeating(temp, cur) == False : + print "Transition to cooling (i=%d)" % i + heating = False + i += 1 + elif heating == False and isHeating(temp, cur) == True : + print "Transition to heating (i=%d)" % i + heating = True + i += 1 + if verbose : + print " Restoring gains" + self.setHeatingGains(*orig_heat_gains) + self.setCoolingGains(*orig_cool_gains) + if log == True : + log = 1# MARK + +def _test_tempController() : + t = tempController(controller=1, maxCurrent=0.1) + + print "Temp = %g" % t.getTemp() + print "Current = %g" % t.getCurrent() + print "Setpoint = %g" % t.getSetpoint() + + print "Setting setpoint to 5.0 deg C" + t.setSetpoint(5.0) + sp = t.getSetpoint() + print "Setpoint = %g" % sp + if sp != 5.0 : + raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0) + time.sleep(10) # give the controller some time to overcome any integral gain + c = t.getCurrent() + print "Current = %g" % c + mca, mcb, mct = t.getMaxCurrent() + if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb + if c != mcb : + raise Exception, "Current not at max %g, and we're shooting for a big temp" % mcb + else : + if c != mca : + raise Exception, "Current not at max %g, and we're shooting for a big temp" % mca + + + print "Setting setpoint to 50.0 deg C" + t.setSetpoint(50.0) + sp = t.getSetpoint() + print "Setpoint = %g" % sp + if sp != 5.0 : + raise Exception, "Setpoint in %g != setpoint out %g" % (sp, 5.0) + time.sleep(10) + c = t.getCurrent() + print "Current = %g" % c + print "Success" + mca, mcb, mct = t.getMaxCurrent() + if t.getTemp() < mct : # we're below the high power limit setpoint, use mcb + if -c != mcb : + raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mcb) + else : + if -c != mca : + raise Exception, "Current not at min %g, and we're shooting for a big temp" % (-mca) + +def test() : + _test_tempController() + +if __name__ == "__main__" : + test() -- 2.26.2