- refactoring of the closedloopunitidentifier

Dynamic/Identification/ClosedLoopGainGlobalSearchResults.cs

@@ -25,6 +25,11 @@ internal class ClosedLoopGainGlobalSearchResults
         /// list of covariance between d_Est and y_Set, calculated for each linear gains
         /// </summary>
         public List<double> covBtwDestAndYsetList;
+        /// <summary>
+        /// self-variance of d_est, calculated for each linear gains
+        /// </summary>
+        public List<double> uPidVarianceList;
+
         /// <summary>
         /// self-variance of d_est, calculated for each linear gains
         /// </summary>
@@ -47,21 +52,23 @@ internal class ClosedLoopGainGlobalSearchResults
         internal ClosedLoopGainGlobalSearchResults()
         {
             unitParametersList= new List<UnitParameters>();
+            uPidVarianceList = new List<double>();
             dEstVarianceList = new List<double>();
             covBtwDestAndYsetList = new List<double>();
             pidLinearProcessGainList = new List<double>();
             covBtwDestAndUexternal = new List<double>();
         }
 
-        public void Add(double gain, UnitParameters unitParameters, double covBtwDestAndYset, double dest_variance,
+        public void Add(double gain, UnitParameters unitParameters, double covBtwDestAndYset, double upidVariance,double eEstVariance,
             double covBtwDestAndUexternal)
         {
           //  var newParams = unitParameters.CreateCopy();
 
             pidLinearProcessGainList.Add(gain);
             unitParametersList.Add(unitParameters);
             covBtwDestAndYsetList.Add(covBtwDestAndYset);
-            dEstVarianceList.Add(dest_variance);
+            uPidVarianceList.Add(upidVariance);
+            dEstVarianceList.Add(eEstVariance);
             this.covBtwDestAndUexternal.Add(covBtwDestAndUexternal);
         }
 
@@ -123,19 +130,30 @@ double[] Scale(double[] v_in)
             }
             Vec vec = new Vec();
 
-            var v1 = dEstVarianceList.ToArray();
+            var v1 = uPidVarianceList.ToArray();
             var v2 = covBtwDestAndYsetList.ToArray();
             var v3 = covBtwDestAndUexternal.ToArray();
             (double v1_Strength, int min_ind) = MinimumStrength(v1);
             (double v2_Strength, int min_ind_v2) = MinimumStrength(v2);
             (double v3_Strength,int min_ind_v3) = MinimumStrength(v3);
-            // add a scaled v2 to the objective only when v1 is very flat around the minimum
+
+            var v4 = dEstVarianceList.ToArray();
+            (double v4_Strength, int min_ind_v4) = MinimumStrength(v4);
+            // add a scaled v4 to the objective only when v1 is very flat around the minimum
             // as a "tiebreaker"
 
             if (v1_Strength == 0 && v2_Strength == 0 && v3_Strength == 0)
-            {
-                //defaultModel.modelParameters.AddWarning(UnitdentWarnings.ClosedLoopEst_GlobalSearchFailedToFindLocalMinima);
-                return new Tuple<UnitModel, bool>(null, false);
+            {/*
+                if (v4_Strength > 0)
+                {
+                    return  new Tuple<UnitModel, bool>(new UnitModel(unitParametersList.ElementAt(min_ind_v4)), true);
+                }
+                else*/
+                {
+                    return new Tuple<UnitModel, bool>(null, false);
+                }
+                //modelParameters.AddWarning(UnitdentWarnings.ClosedLoopEst_GlobalSearchFailedToFindLocalMinima);
+
             }
             double[] objFun = v1;
 

Dynamic/PlantSimulator/PlantSimulator.cs

@@ -399,15 +399,15 @@ public string ConnectModels(ISimulatableModel upstreamModel, ISimulatableModel d
         /// <summary>
         /// Create a PlantSimulator and TimeSeriesDataSet from a UnitDataSet, PidModel and UnitModel to do closed-loop simulations
         /// <para>
-        /// The feedback loop has no disturbance signal added, but this can be added to the returned PlantSimulator as needed.
+        /// The feedback loop has NO disturbance signal added, but this can be added to the returned PlantSimulator as needed.
         /// </para>
         /// </summary>
         /// <param name="unitDataSet"></param>
         /// <param name="pidModel"></param>
         /// <param name="unitModel"></param>
         /// <param name="pidInputIdx"></param>
-        /// <returns></returns>
-        public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoop(UnitDataSet unitDataSet, PidModel pidModel, 
+        /// <returns>a simulator object and a dataset object that is ready to be simulated with Simulate() </returns>
+        public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoopNoDisturbance(UnitDataSet unitDataSet, PidModel pidModel, 
             UnitModel unitModel, int pidInputIdx=0)
         {
             var plantSim = new PlantSimulator(
@@ -416,8 +416,8 @@ public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoop(UnitDataSet
             var signalId2 = plantSim.ConnectModels(pidModel, unitModel, pidInputIdx);
 
             var inputData = new TimeSeriesDataSet();
-            inputData.Add(signalId1, unitDataSet.Y_meas);
-            inputData.Add(signalId2, unitDataSet.U.GetColumn(pidInputIdx));
+            inputData.Add(signalId1, (double[])unitDataSet.Y_meas.Clone());
+            inputData.Add(signalId2, (double[])unitDataSet.U.GetColumn(pidInputIdx).Clone());
 
             if (unitDataSet.U.GetNColumns() > 1)
             {
@@ -426,17 +426,34 @@ public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoop(UnitDataSet
                     if (curColIdx == pidInputIdx)
                         continue;
                     inputData.Add(plantSim.AddExternalSignal(unitModel, SignalType.External_U, curColIdx),
-                        unitDataSet.U.GetColumn(curColIdx));
+                        (double[])unitDataSet.U.GetColumn(curColIdx).Clone());
                 }
             }
 
-            inputData.Add(plantSim.AddExternalSignal(pidModel, SignalType.Setpoint_Yset), unitDataSet.Y_setpoint);
+            inputData.Add(plantSim.AddExternalSignal(pidModel, SignalType.Setpoint_Yset), (double[])unitDataSet.Y_setpoint.Clone());
             inputData.CreateTimestamps(unitDataSet.GetTimeBase());
             inputData.SetIndicesToIgnore(unitDataSet.IndicesToIgnore);
 
             return (plantSim, inputData);
         }
 
+        /// <summary>
+        /// Create a feedback loop, where the process model has an additive disturbance that is to be estimated.
+        /// </summary>
+        /// <param name="unitDataSet"></param>
+        /// <param name="pidModel"></param>
+        /// <param name="unitModel"></param>
+        /// <param name="pidInputIdx"></param>
+        /// <returns>a simulator object and a dataset object that is ready to be simulated with Simulate() </returns>
+        public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoopWithEstimatedDisturbance(UnitDataSet unitDataSet, PidModel pidModel,
+             UnitModel unitModel, int pidInputIdx = 0)
+        {
+            // vital that signal follows naming convention, otherwise it will not be estimated, but should be provided.
+            unitModel.AddSignalToOutput(SignalNamer.EstDisturbance(unitModel));
+            (PlantSimulator sim, TimeSeriesDataSet data) =  CreateFeedbackLoopNoDisturbance(unitDataSet,pidModel, unitModel, pidInputIdx);
+            return (sim,data);
+        }
+
 
         /// <summary>
         /// Get a TimeSeriesDataSet of all external signals of model.
@@ -593,6 +610,7 @@ static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl
         ///  If the model is a unitModel and the inputData inludes both the measured y and measured u, the
         ///  simData will include an estimate of the additive disturbance.
         ///  </para>
+        /// 
         /// </summary>
         /// <param name="inputData"></param>
         /// <param name="singleModelName"></param>
@@ -707,7 +725,15 @@ public bool SimulateSingle(TimeSeriesDataSet inputData, string singleModelName,
             return true;
         }
         /// <summary>
-        /// Perform a dynamic simulation of the model provided, given the specified connections and external signals. 
+        /// Perform a "plant-wide" full dynamic simulation of the entire plant,i.e. all models in the plant, given the specified connections and external signals. 
+        /// <para>
+        /// The dynamic simulation will also return estimated disturbances in simData, if the plant contains feedback loops where there is an additive 
+        /// disturbance with a signal named according to SignalNamer.EstDisturbance() convention
+        /// </para>
+        /// <para>
+        ///  The simulation will also set the <c>PlantFitScore</c> which can be used to evalute the fit of the simulation to the plant data.
+        ///  For this score to be calculated, the measured time-series corresponding to <c>simData</c> need to be provided in <c>inputData</c>
+        ///  </para>
         /// </summary>
         /// <param name="inputData">the external signals for the simulation(also, determines the simulation time span and timebase)</param>
         /// <param name="simData">the simulated data set to be outputted(excluding the external signals)</param>

Dynamic/PlantSimulator/SignalNamer.cs

@@ -22,6 +22,11 @@ public class SignalNamer
         /// <returns>a unique string identifier that is used to identify a signal</returns>
         public static string GetSignalName(string modelID, SignalType signalType, int idx = 0)
         {
+            if (signalType == SignalType.Disturbance_D)
+            {
+                return EstDisturbance(modelID);
+            }
+
             if (idx == 0)
                 return modelID + separator + signalType.ToString();
             else

Dynamic/SimulatableModels/UnitModel.cs

@@ -67,6 +67,9 @@ public class UnitModel : ModelBaseClass, ISimulatableModel
 
         private List<ProcessTimeDelayIdentWarnings> TimeDelayEstWarnings { get; }
 
+        /// <summary>
+        /// Empty constructor
+        /// </summary>
         public UnitModel(){}
 
         /// <summary>
@@ -93,10 +96,7 @@ public UnitModel(UnitParameters modelParameters, UnitDataSet dataSet,string ID =
             this.ID = ID;
             InitSim(modelParameters);
         }
-        public string test()
-        {
-            return "test";
-        }
+
         /// <summary>
         /// Answers if the model can be simulated with the inputs provided.
         /// </summary>

TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs

@@ -112,7 +112,7 @@ public void PlantSimulator_StepDisturbance_EstimatesOk(double stepAmplitude)
         [TestCase(4,-1)]
         public void PlantSimulator_StepDisturbanceANDSetPointStep_EstimatesOk(double disturbanceStepAmplitude,double setpointAmplitude)
         {
-           // Shared.EnablePlots();
+          //  Shared.EnablePlots();
             var trueDisturbance = TimeSeriesCreator.Step(100, N, 0, disturbanceStepAmplitude);
             var setpoint = TimeSeriesCreator.Step(50, N, 50, 50+setpointAmplitude);
             DisturbanceTestUsingPlantSimulator(new UnitModel(dynamicModelParameters, "PlantSim_d"), trueDisturbance, true, setpoint);
@@ -181,9 +181,17 @@ public void DisturbanceTestUsingPlantSimulator(UnitModel processModel, double[]
             // 2.create plant model without disturbance, and try to to find the disturbance signal
             {
                 var pidModel1 = new PidModel(pidParameters1, "PID1");
-              
+
                 // important to add a disturbance signal with naming convention to the output of the process, to signal that the disturbance is to be estimated!
-                   var processModel3 = new UnitModel(dynamicModelParameters, "Proc1");
+
+                // begin testing
+                // var dynamicModelParameters_new = dynamicModelParameters.CreateCopy();
+                // dynamicModelParameters_new.TimeConstant_s = 20;
+                //end testing
+
+                //var processModel3 = new UnitModel(dynamicModelParameters_new, "Proc1");
+                var processModel3 = new UnitModel(dynamicModelParameters, "Proc1");
+
                 var distSignal = SignalNamer.EstDisturbance(processModel3);
                 processModel3.AddSignalToOutput(distSignal);
                 var plantSim = new PlantSimulator(
@@ -214,9 +222,10 @@ public void DisturbanceTestUsingPlantSimulator(UnitModel processModel, double[]
                 inputData.CreateTimestamps(timeBase_s);
 
                 //////////////////////////////////
-                var isOK = plantSim.Simulate(inputData, out TimeSeriesDataSet simDataSetWithDisturbance);
+                var isOK = plantSim.Simulate(inputData, out var simDataSetWithDisturbance);
                 //////////////////////////////////
                 Assert.IsTrue(isOK);
+                Assert.IsTrue(plantSim.PlantFitScore == 100);
                 Assert.IsTrue(simDataSetWithDisturbance.ContainsSignal(distSignal));
                 if (doAssertResult)
                 {

TimeSeriesAnalysis.Tests/Tests/FindDisturbanceAndModelSimultanouslyTester_SISO.cs

@@ -27,7 +27,7 @@ class FindDisturbanceAndModelSimultanouslyTester_SISO
         {
             TimeConstant_s = 10,
             LinearGains = new double[] { 1.5 },
-            TimeDelay_s = 5,
+            TimeDelay_s = 0,//was: 5
             Bias = 5
         };
 
@@ -170,9 +170,16 @@ public void Dynamic_SetpointStep( double ysetStepAmplitude)
                 false, true, yset, precisionPrc);
         }
 
+        /*
+        I think the reason this test struggles is that closedloopestimator tries to find the process model that results in the
+        disturbance with the smallest average change, but for continously acting disturbances like this sinus disturbance, 
+        this may not be as good an assumption as for the step disturbances considered in other tests. 
+        */
+
         [TestCase(5, 1.0, Category="NotWorking_AcceptanceTest"), NonParallelizable]
         [TestCase(1, 1.0, Category = "NotWorking_AcceptanceTest") ]
-        [TestCase(1, 5.0)]
+        [TestCase(1, 5.0)]// this only works when the step change is much bigger than the disturbance
+        [TestCase(1, 0.0, Category = "NotWorking_AcceptanceTest")]
         public void Static_SinusDistANDSetpointStep(double distSinusAmplitude,
             double ysetStepAmplitude)
         {
@@ -257,9 +264,9 @@ public void StepAtStartOfDataset_IsExcludedFromAnalysis()
                 trueDisturbance, false, true,null,10, doAddBadData);
         }
 
-        [TestCase(-5,5)]
-        [TestCase(5, 5)]
-        [TestCase(10, 5)]
+        [TestCase(-5,10)]
+        [TestCase(5, 10)]
+        [TestCase(10, 10)]
         public void Dynamic_DistStep_EstimatesOk(double stepAmplitude, double processGainAllowedOffsetPrc, 
             bool doNegativeGain =false)
         {
@@ -311,10 +318,13 @@ public void GenericDisturbanceTest  (UnitModel trueProcessModel, double[] trueDi
                 pidDataSet.Y_setpoint[50] = Double.NaN;
                 pidDataSet.Y_meas[400] = Double.NaN;
                 pidDataSet.U[500,0] = Double.NaN;
-
             }
+            // NB! uses the "perfect" pid-model in the identification process
+
+            var estPidParam = new PidParameters(pidModel1.GetModelParameters());
+          //  estPidParam.Kp = estPidParam.Kp * 0.5;// try adding a wrong kp and see what it does
 
-            (var identifiedModel, var estDisturbance) = ClosedLoopUnitIdentifier.Identify(pidDataSet, pidModel1.GetModelParameters());
+            (var identifiedModel, var estDisturbance) = ClosedLoopUnitIdentifier.Identify(pidDataSet, estPidParam);
 
             Console.WriteLine(identifiedModel.ToString());
             Console.WriteLine();

TimeSeriesAnalysis/TimeSeriesDataSet.cs

@@ -536,7 +536,7 @@ public double[] GetValues(string signalName)
                 return Vec<double>.Fill(dataset_constants[signalName],N.Value);   
             else
             {
-                Shared.GetParserObj().AddError("TimeSeriesDataSert.GetValues() did not find signal:" + signalName);
+                Shared.GetParserObj().AddError("TimeSeriesDataSet.GetValues() did not find signal:" + signalName);
                 return null;
             }
         }

Utility/PlotGain.cs

@@ -183,6 +183,7 @@ public static void PlotGainSched(GainSchedModel model1, GainSchedModel model2 =
         /// </summary>
         /// <param name="model"></param>
         /// <param name="modelName"></param>
+        /// <param name="comment"></param>
         /// <param name="inputIdx"></param>
         /// <param name="uMinExt"></param>
         /// <param name="uMaxExt"></param>