- working on refactoring the DisturbanceEstimator

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
         /// <param name="pidParams">an estimate of the PID-parameters</param>
         /// <param name="pidInputIdx">the index of U in the above dataset that is driven by the PID-controller.</param>
         /// <returns></returns>
-        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<ISimulatableModel> { 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<double[]> { unitData.Y_sim,unitData.Y_meas, 
-                    unitData.U_sim.GetColumn(0),  unitData.U.GetColumn(0)},
-                    new List<string> { "y1=y_sim", "y1=y_meas", "y3=u_sim","y3=u_meas" },
-                    unitData.GetTimeBase(), "doDebuggingPlot" + "_closedlooopsim_"+name); 
-            }
-
-            return isOk;*/
         }
     }
 }

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<double[]> { d_est, d_est2, d_HF }, new List<string> { "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;

TimeSeriesAnalysis.Tests/Tests/CluiCommonTests.cs

@@ -103,23 +103,28 @@ 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();
 
             double distTrueAmplitude = vec.Max(vec.Abs(trueDisturbance));
             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<double[]>{ pidDataSet.Y_meas, pidDataSet.Y_setpoint,pidDataSet.Y_proc,
-                pidDataSet.U.GetColumn(0),  estDisturbance, trueDisturbance },
-                    new List<string> { "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<double[]>{ 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<string> { "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();
             }
 

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);
 
 

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})