From b81d9686cc8f981461666824fca8291ad9170704 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steinar=20Elgs=C3=A6ter?= Date: Thu, 12 Dec 2024 14:57:51 +0100 Subject: [PATCH] - working on refactoring the DisturbanceEstimator --- .../ClosedLoopUnitIdentifier.cs | 29 ++----------------- .../Identification/DisturbanceIdentifier.cs | 18 +++++++++--- .../Tests/CluiCommonTests.cs | 19 +++++++----- .../Tests/DisturbanceEstimatorTests.cs | 2 +- docs/articles/sysid_disturbance.md | 19 ++++++++++++ 5 files changed, 48 insertions(+), 39 deletions(-) diff --git a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs index c74fbab..2dad165 100644 --- a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs +++ b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs @@ -120,7 +120,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters // run1: no process model assumed, let disturbance estimator guesstimate a pid-process gain, // to give afirst estimate of the disturbance - var unitModel_step1 = EstimateClosedLoopProcessGain(dataSetRun1, pidParams, pidInputIdx); + var unitModel_step1 = ModelFreeEstimateClosedLoopProcessGain(dataSetRun1, pidParams, pidInputIdx); var distIdResult1 = DisturbanceIdentifier.CalculateDisturbanceVector(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); if (doConsoleDebugOut) Console.WriteLine("Step1: " + unitModel_step1.GetModelParameters().LinearGains.ElementAt(pidInputIdx).ToString("F3", CultureInfo.InvariantCulture)); @@ -901,11 +901,8 @@ private static double GuessSignOfProcessGain(UnitDataSet unitDataSet, PidParamet /// an estimate of the PID-parameters /// the index of U in the above dataset that is driven by the PID-controller. /// - private static UnitModel EstimateClosedLoopProcessGain(UnitDataSet unitDataSet, PidParameters pidParams, int pidInputIdx) + private static UnitModel ModelFreeEstimateClosedLoopProcessGain(UnitDataSet unitDataSet, PidParameters pidParams, int pidInputIdx) { - - - var unitModel = new UnitModel(); var vec = new Vec(unitDataSet.BadDataID); @@ -1066,29 +1063,7 @@ public static bool ClosedLoopSim(UnitDataSet unitData, UnitParameters modelParam unitData.Y_sim = simData.GetValues(unitModel.GetID(), SignalType.Output_Y, 0); unitData.Y_proc = simData.GetValues(unitModel.GetID()); } - return isOk; - - /* - - - - // TODO: replace this with a "closed-loop" simulator that uses the PlantSimulator. - //var ps = new PlantSimulator(new List { unitModel,pidModel }); - //var inputData = new TimeSeriesDataSet();// TODO: need to map signal names here. - //ps.Simulate(inputData, out var simData); - - var sim = new UnitSimulator(unitModel); - bool isOk = sim.CoSimulate(pidModel, ref unitData); - if (doDebuggingPlot) - { - Plot.FromList(new List { unitData.Y_sim,unitData.Y_meas, - unitData.U_sim.GetColumn(0), unitData.U.GetColumn(0)}, - new List { "y1=y_sim", "y1=y_meas", "y3=u_sim","y3=u_meas" }, - unitData.GetTimeBase(), "doDebuggingPlot" + "_closedlooopsim_"+name); - } - - return isOk;*/ } } } \ No newline at end of file diff --git a/Dynamic/Identification/DisturbanceIdentifier.cs b/Dynamic/Identification/DisturbanceIdentifier.cs index 90a8a51..c56627d 100644 --- a/Dynamic/Identification/DisturbanceIdentifier.cs +++ b/Dynamic/Identification/DisturbanceIdentifier.cs @@ -293,9 +293,9 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat // non-disturbance related changes in the dataset producing "unitDataSet_adjusted" var unitDataSet_adjusted = RemoveSetpointAndOtherInputChangeEffectsFromDataSet(unitDataSet, unitModel, pidInputIdx, pidParams); unitDataSet_adjusted.D = null; - (bool isOk, double[] y_sim) = PlantSimulator.SimulateSingle(unitDataSet_adjusted, unitModel, false); + (bool isOk, double[] y_proc) = PlantSimulator.SimulateSingle(unitDataSet_adjusted, unitModel, false); - if (y_sim == null) + if (y_proc == null) { result.zeroReason = DisturbanceSetToZeroReason.UnitSimulatorUnableToRun; return result; @@ -314,7 +314,7 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat } } - double[] d_LF = vec.Multiply(vec.Subtract(y_sim, y_sim[indexOfFirstGoodValue]), -1); + double[] d_LF = vec.Multiply(vec.Subtract(y_proc, y_proc[indexOfFirstGoodValue]), -1); double[] d_HF = vec.Subtract(unitDataSet_adjusted.Y_meas, unitDataSet_adjusted.Y_setpoint); // d_u : (low-pass) back-estimation of disturbances by the effect that they have on u as the pid-controller integrates to @@ -324,7 +324,17 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat // will influence the process model identification afterwards. // d = d_HF+d_LF - double[] d_est = vec.Add(d_HF, d_LF); + // double[] d_est2 = vec.Add(d_HF, d_LF); + // alternative: + double[] d_est = vec.Subtract(unitDataSet_adjusted.Y_meas, y_proc); + +// Shared.EnablePlots(); + // Plot.FromList(new List { d_est, d_est2, d_HF }, new List { "y1=d1", "y1=d2","y1=d_HF" }, unitDataSet.Times, "testplot"); + // Shared.DisablePlots(); + + + + result.d_est = d_est; result.d_LF = d_LF; result.d_HF = d_HF; diff --git a/TimeSeriesAnalysis.Tests/Tests/CluiCommonTests.cs b/TimeSeriesAnalysis.Tests/Tests/CluiCommonTests.cs index 8eb1b7d..c8864b9 100644 --- a/TimeSeriesAnalysis.Tests/Tests/CluiCommonTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/CluiCommonTests.cs @@ -103,8 +103,9 @@ static void DisturbancesToString(double[] estDisturbance, double[] trueDisturban - static public void CommonPlotAndAsserts(UnitDataSet pidDataSet, double[] estDisturbance, double[] trueDisturbance, - UnitModel identifiedModel, UnitModel trueModel, double maxAllowedGainOffsetPrc, double maxAllowedMeanDisturbanceOffsetPrc = 30, bool isStatic=false) + static public void CommonPlotAndAsserts(UnitDataSet unitDataSet, double[] estDisturbance, double[] trueDisturbance, + UnitModel identifiedModel, UnitModel trueModel, double maxAllowedGainOffsetPrc, + double maxAllowedMeanDisturbanceOffsetPrc = 30, bool isStatic=false) { Vec vec = new Vec(); @@ -112,14 +113,18 @@ static public void CommonPlotAndAsserts(UnitDataSet pidDataSet, double[] estDist Assert.IsTrue(estDisturbance != null); string caseId = TestContext.CurrentContext.Test.Name.Replace("(", "_"). Replace(")", "_").Replace(",", "_") + "y"; - bool doDebugPlot = true; + bool doDebugPlot = false; if (doDebugPlot) { + double[] d_HF = vec.Subtract(unitDataSet.Y_meas, unitDataSet.Y_setpoint); + double[] d_LF = vec.Multiply(vec.Subtract(unitDataSet.Y_proc, unitDataSet.Y_proc[0]), -1); + Shared.EnablePlots(); - Plot.FromList(new List{ pidDataSet.Y_meas, pidDataSet.Y_setpoint,pidDataSet.Y_proc, - pidDataSet.U.GetColumn(0), estDisturbance, trueDisturbance }, - new List { "y1=y_meas", "y1=y_set", "y1=y_sim", "y2=u(right)", "y3=est disturbance", "y3=true disturbance" }, - pidDataSet.GetTimeBase(), caseId + "commonplotandasserts"); + Plot.FromList(new List{ unitDataSet.Y_meas, unitDataSet.Y_setpoint,unitDataSet.Y_proc, + unitDataSet.U.GetColumn(0), unitDataSet.U_sim.GetColumn(0), estDisturbance, trueDisturbance, d_HF,d_LF }, + new List { "y1=y_meas", "y1=y_set", "y1=y_proc", "y2=u_meas(right)","y2=u_sim(right)", + "y3=est disturbance", "y3=true disturbance","y3= d_HF","y3=d_LF"}, + unitDataSet.GetTimeBase(), caseId + "commonplotandasserts"); Shared.DisablePlots(); } diff --git a/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs b/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs index c044e64..c6b5b03 100644 --- a/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs @@ -41,7 +41,7 @@ class FindDisturbanceWhenProcessModelIsGiven }; int timeBase_s = 1; - int N = 300;// TODO:influences the results! + int N = 300; DateTime t0 = new DateTime(2010,1,1); diff --git a/docs/articles/sysid_disturbance.md b/docs/articles/sysid_disturbance.md index 6b42259..41ae8cc 100644 --- a/docs/articles/sysid_disturbance.md +++ b/docs/articles/sysid_disturbance.md @@ -211,7 +211,26 @@ Some obervations: It could be an idea to count the number of zero crossing of e in the dataset, if it is very low, that may be an indication that the process gain is low. +*It could be important to look at the internal output of the process model.* +*It could be that it is possible to have a very good estiamte of the distrubance d, yet the process gain can still be of by 40 percent. +This has been observed in some random walk unit tests where the process gain is small. This could be because the d_est in this cae is mostly +determined by ``d_HF```? This could be the clue to determining the uncertainty of the process gain, looking at how much the disturbance vector +estimate is sensitive to the process gain values?* + +There are essentially two ways of calculating the disturbance + +1. By subtracting the modelled ``y_{proc}(\hat{u})`` from ``\bar{y}'' : ``d_est = \bar{y} - y_{proc}(\hat{u})`` +2. By `d_est = d_HF(\hat{u}, y_{set}) +d_LF (u)`` + +where +``d_LF (u) = \hat{y}(u(t))- \hat{y} (u(t_0))`` + +Note that ``d_HF`` does not change with changing estimtes of the model gain or other paramters, while d_LF does. + +Combining the two above means that + +``d_LF(u) = (\bar{y} - y_{proc}(\hat{u})) - d_HF(\hat{u}, y_{set})