Version 1.15rc3

Improved PID tuning

7 files changed, 44 insertions, 43 deletions

diff --git a/src/Configuration.h b/src/Configuration.h
index c75a98b6..def6e1f9 100644
--- a/src/Configuration.h
+++ b/src/Configuration.h
@@ -26,7 +26,7 @@ Licence: GPL
 // Firmware name is now defined in the Pins file
 
 #ifndef VERSION
-# define VERSION "1.15-rc2"
+# define VERSION "1.15-rc3"
 #endif
 
 #ifndef DATE
diff --git a/src/GCodes/GCodes.cpp b/src/GCodes/GCodes.cpp
index 93885324..371fe517 100644
--- a/src/GCodes/GCodes.cpp
+++ b/src/GCodes/GCodes.cpp
@@ -3206,7 +3206,7 @@ bool GCodes::HandleMcode(GCodeBuffer* gb, StringRef& reply)
 				if (!success)
 				{
 					error = true;
-					platform->MessageF(GENERIC_MESSAGE, "Setting pin %d to %d is not supported\n", pin, val);
+					reply.printf("Setting pin %d to %d is not supported\n", pin, val);
 				}
 			}
 		}
@@ -4216,9 +4216,13 @@ bool GCodes::HandleMcode(GCodeBuffer* gb, StringRef& reply)
 					{
 						// When reporting the PID parameters, we scale them by 255 for compatibility with older firmware and other firmware
 						const PidParams& spParams = model.GetPidParameters(false);
-						reply.catf("\nSetpoint change: P%.1f, I%.2f, D%.1f", 255.0 * spParams.kP, 255.0 * spParams.kI, 255.0 * spParams.kD);
+						const float scaledSpKp = 255.0 * spParams.kP;
+						reply.catf("\nSetpoint change: P%.1f, I%.2f, D%.1f",
+								scaledSpKp, scaledSpKp * spParams.recipTi, scaledSpKp * spParams.tD);
 						const PidParams& ldParams = model.GetPidParameters(true);
-						reply.catf("\nLoad change: P%.1f, I%.2f, D%.1f", 255.0 * ldParams.kP, 255.0 * ldParams.kI, 255.0 * ldParams.kD);
+						const float scaledLoadKp = 255.0 * ldParams.kP;
+						reply.catf("\nLoad change: P%.1f, I%.2f, D%.1f",
+								scaledLoadKp, scaledLoadKp * ldParams.recipTi, scaledLoadKp * ldParams.tD);
 					}
 				}
 			}
diff --git a/src/Heating/FOPDT.cpp b/src/Heating/FOPDT.cpp
index 2682c20f..6515cb12 100644
--- a/src/Heating/FOPDT.cpp
+++ b/src/Heating/FOPDT.cpp
@@ -57,26 +57,13 @@ bool FopDt::SetParameters(float pg, float ptc, float pdt, float pMaxPwm, bool pU
 void FopDt::CalcPidConstants()
 {
 	const float timeFrac = deadTime/timeConstant;
-
-#if 1
 	loadChangeParams.kP = 0.7/(gain * timeFrac);
-	loadChangeParams.kI = loadChangeParams.kP/(deadTime * 2.0);
-	loadChangeParams.kD = loadChangeParams.kP * deadTime;
+	loadChangeParams.recipTi = 1.0/(deadTime * 2.0);								// Ti = 2 * deadTime
+	loadChangeParams.tD = deadTime * 0.7;
 
 	setpointChangeParams.kP = 0.7/(gain * timeFrac);
-	setpointChangeParams.kI = setpointChangeParams.kP/max<float>(deadTime * 2.0, timeConstant);
-	setpointChangeParams.kD = setpointChangeParams.kP * deadTime;
-#else
-	// Calculate the PID parameters for responding to changes in load, using ITAE-load
-	loadChangeParams.kP = (0.77902/gain) * pow(timeFrac, -1.06401);
-	loadChangeParams.kI = loadChangeParams.kP/((1.0/1.14311) * timeConstant * pow(timeFrac, 0.70949));
-	loadChangeParams.kD = loadChangeParams.kP * 0.57137 * timeConstant * pow(timeFrac, 1.03826);
-
-	// Calculate the PID parameters for responding to changes in the setpoint using IAE-setpoint
-	setpointChangeParams.kP = (0.65/gain) * pow(timeFrac, -1.04432);
-	setpointChangeParams.kI = setpointChangeParams.kP * (0.9895 + 0.09539 * timeFrac)/timeConstant;
-	setpointChangeParams.kD = setpointChangeParams.kP * 0.50814 * timeConstant * pow(timeFrac, 1.08433);
-#endif
+	setpointChangeParams.recipTi = 1.0/max<float>(deadTime * 2.0, timeConstant);	// Ti = time constant unless dead time is very long
+	setpointChangeParams.tD = deadTime * 0.7;
 }
 
 // End
diff --git a/src/Heating/FOPDT.h b/src/Heating/FOPDT.h
index 82688076..d174da53 100644
--- a/src/Heating/FOPDT.h
+++ b/src/Heating/FOPDT.h
@@ -12,9 +12,9 @@
 
 struct PidParams
 {
-	float kP;
-	float kI;
-	float kD;
+	float kP;			// controller (not model) gain
+	float recipTi;		// reciprocal of controller integral time
+	float tD;			// controller differential time
 };
 
 class FopDt
diff --git a/src/Heating/Heat.cpp b/src/Heating/Heat.cpp
index 53fb0220..49ca6c67 100644
--- a/src/Heating/Heat.cpp
+++ b/src/Heating/Heat.cpp
@@ -134,8 +134,9 @@ Heat::HeaterStatus Heat::GetStatus(int8_t heater) const
 
 	return (pids[heater]->FaultOccurred()) ? HS_fault
 			: (pids[heater]->SwitchedOff()) ? HS_off
-				: (pids[heater]->Active()) ? HS_active
-					: HS_standby;
+				: (pids[heater]->IsTuning()) ? HS_tuning
+					: (pids[heater]->Active()) ? HS_active
+						: HS_standby;
 }
 
 void Heat::SetActiveTemperature(int8_t heater, float t)
diff --git a/src/Heating/Heat.h b/src/Heating/Heat.h
index 00fb9ed5..49dff9df 100644
--- a/src/Heating/Heat.h
+++ b/src/Heating/Heat.h
@@ -31,7 +31,7 @@ class Heat
 {
 public:
 	// Enumeration to describe the status of a heater. Note that the web interface returns the numerical values, so don't change them.
-	enum HeaterStatus { HS_off = 0, HS_standby = 1, HS_active = 2, HS_fault = 3 };
+	enum HeaterStatus { HS_off = 0, HS_standby = 1, HS_active = 2, HS_fault = 3, HS_tuning = 4 };
 
 	Heat(Platform* p);
 	void Spin();												// Called in a tight loop to keep everything going
diff --git a/src/Heating/Pid.cpp b/src/Heating/Pid.cpp
index 2b49dcf0..bf13b219 100644
--- a/src/Heating/Pid.cpp
+++ b/src/Heating/Pid.cpp
@@ -257,26 +257,30 @@ void PID::Spin()
 			if (usingPid)
 			{
 				// Using PID mode
-				float kP, kI, kD, gain;
+				float kP, recipTi, tD, gain;
+				bool inSetPointMode;
 				if (useModel)
 				{
-					const PidParams& params = model.GetPidParameters(fabs(error) < 1.0);
+					inSetPointMode = fabs(error) > 1.0;			// use modified PID when the error is large
+					const PidParams& params = model.GetPidParameters(!inSetPointMode);
 					kP = params.kP;
-					kI = params.kI;
-					kD = params.kD;
+					recipTi = params.recipTi;
+					tD = params.tD;
 					gain = model.GetGain();
 				}
 				else
 				{
-					kP = pp.kP/255.0 * pp.kS;
-					kI = pp.kI/255.0 * pp.kS;
-					kD = pp.kD/255.0 * pp.kS;
+					inSetPointMode = false;						// use standard PID always
+					kP = (pp.kP * pp.kS) * (1.0/255.0);
+					recipTi = pp.kI/pp.kP;
+					tD = pp.kD/kP;
 					gain = 255.0/pp.kT;
 				}
 
 				// If the P term still has full authority, preset the I term to the expected PWM for this temperature
 				// and turn the heater full on or full off
-				const float pPlusD = (kP * error) - (kD * derivative);
+				const float errorMinusDterm = error - (tD * derivative);
+				const float pPlusD = kP * errorMinusDterm;
 				const float expectedPwm = constrain<float>((temperature - NormalAmbientTemperature)/gain, 0.0, maxPwm);
 				if (pPlusD + expectedPwm > maxPwm)
 				{
@@ -290,7 +294,12 @@ void PID::Spin()
 				}
 				else
 				{
-					iAccumulator = constrain<float>(iAccumulator + (error * kI * platform->HeatSampleInterval() * MillisToSeconds), 0.0, maxPwm);
+					// In the following we use a modified PID when the temperature is a long way off target.
+					// During initial heating or cooling, the D term represents expected overshoot, which we don't want to add to the I accumulator.
+					// When we are in load mode, the I term is much larger and the D term doesn't represent overshoot, so use normal PID.
+					const float errorToUse = (inSetPointMode) ? errorMinusDterm : error;
+					iAccumulator = constrain<float>(iAccumulator + (errorToUse * kP * recipTi * platform->HeatSampleInterval() * MillisToSeconds),
+													0.0, maxPwm);
 					lastPwm = constrain<float>(pPlusD + iAccumulator, 0.0, maxPwm);
 				}
 			}
@@ -549,7 +558,6 @@ void PID::DoTuningStep()
 				timeSetHeating = tuningPhaseStartTime = millis();
 				lastPwm = tuningPwm;										// turn on heater at specified power
 				tuningReadingInterval = platform->HeatSampleInterval();		// reset sampling interval
-				active = true;
 				mode = HeaterMode::tuning1;
 				return;
 			}
@@ -580,10 +588,10 @@ void PID::DoTuningStep()
 						DisplayBuffer("At phase 1 end");
 					}
 					tuningTimeOfFastestRate = index * tuningReadingInterval * MillisToSeconds;
-					tuningPhaseStartTime += index * tuningReadingInterval;
-					tuningFastestRate = (tuningTempReadings[index + 1] - tuningTempReadings[index - 1]) / (tuningReadingInterval * 2 * MillisToSeconds);
+					tuningFastestRate = (tuningTempReadings[index + 2] - tuningTempReadings[index - 2]) / (tuningReadingInterval * 4 * MillisToSeconds);
 
 					// Move the readings down so as to start at the max rate index
+					tuningPhaseStartTime += index * tuningReadingInterval;
 					tuningReadingsTaken -= index;
 					for (size_t i = 0; i < tuningReadingsTaken; ++i)
 					{
@@ -599,7 +607,7 @@ void PID::DoTuningStep()
 			{
 				// In the following, the figure of 2.75 was chosen because a value of 2 is too low to handle the bed heater
 				// with embedded thermistor on my Kossel (reservoir effect)
-				if (ReadingsStable(tuningReadingsTaken/2, (temperature - tuningTempReadings[0]) * 0.269))		// if we have been going for ~2 time constants
+				if (ReadingsStable(tuningReadingsTaken/2, (temperature - tuningTempReadings[0]) * 0.2))		// if we have been going for ~2.75 time constants
 				{
 					// We have been heating for long enough, so we can do a fit
 					FitCurve();
@@ -643,13 +651,14 @@ bool PID::ReadingsStable(size_t numReadings, float maxDiff) const
 }
 
 // Return the index in the temperature readings of the maximum rate of increase
+// In view of the poor resolution of most thermistors at low temperatures, we measure over 4 time intervals instead of 2.
 /*static*/ size_t PID::GetMaxRateIndex()
 {
-	size_t index = 1;
+	size_t index = 2;
 	float maxIncrease = 0.0;
-	for (size_t i = 1; i + 1 < tuningReadingsTaken; ++i)
+	for (size_t i = 2; i + 2 < tuningReadingsTaken; ++i)
 	{
-		const float increase = tuningTempReadings[i + 1] - tuningTempReadings[i - 1];
+		const float increase = tuningTempReadings[i + 2] - tuningTempReadings[i - 2];
 		if (increase > maxIncrease)
 		{
 			maxIncrease = increase;