From 6516a666332763d0dcf72db7347a9a813fd3fc6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steinar=20Elgs=C3=A6ter?= Date: Tue, 17 Dec 2024 10:17:50 +0100 Subject: [PATCH] -renaming disturbanceidentifier->disturbancecalculator -fix issues with PlantSimulator.SimulateSingle when using it to back-calculate the disturbance --- .../ClosedLoopUnitIdentifier.cs | 15 +++--- ...Identifier.cs => DisturbanceCalculator.cs} | 14 ++--- Dynamic/PlantSimulator/PlantSimulator.cs | 52 +++++++++++++++---- .../Tests/DisturbanceEstimatorTests.cs | 5 +- docs/articles/processsimulator.md | 2 +- 5 files changed, 60 insertions(+), 28 deletions(-) rename Dynamic/Identification/{DisturbanceIdentifier.cs => DisturbanceCalculator.cs} (97%) diff --git a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs index bb1c4d2..0100c46 100644 --- a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs +++ b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs @@ -118,12 +118,13 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters var GSdescription = ""; - // ---------------- // step 1: no process model assumed, let disturbance estimator guesstimate a pid-process gain, // to give a first estimate of the disturbance var unitModel_step1 = ModelFreeEstimateClosedLoopProcessGain(dataSetRun1, pidParams, pidInputIdx); - var distIdResult1 = DisturbanceIdentifier.CalculateDisturbanceVector(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); + var distIdResult1 = DisturbanceCalculator.CalculateDisturbanceVector(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); + if(distIdResult1.ErrorType != DisturbanceEstimationError.None) + Console.WriteLine("Step 1 Error:" + distIdResult1.ErrorType); if (doConsoleDebugOut) ConsoleDebugOut(unitModel_step1, "Step 1"); idUnitModelsList.Add(unitModel_step1); @@ -312,7 +313,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters // approach1: // "look for the time-constant that gives the disturbance that changes the least" // - var distIdResult2 = DisturbanceIdentifier.CalculateDisturbanceVector + var distIdResult2 = DisturbanceCalculator.CalculateDisturbanceVector (dataSetStep3, unitModel, pidInputIdx, pidParams); var estDisturbances = new List(); var distDevs = new List(); @@ -329,7 +330,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters var newParams = unitModel.GetModelParameters().CreateCopy(); newParams.TimeConstant_s = candiateTc_s; var newModel = new UnitModel(newParams); - var distIdResult_Test = DisturbanceIdentifier.CalculateDisturbanceVector + var distIdResult_Test = DisturbanceCalculator.CalculateDisturbanceVector (dataSetStep3, newModel, pidInputIdx, pidParams); estDisturbances.Add(distIdResult_Test.d_est); @@ -357,7 +358,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters step3params.TimeConstant_s = candiateTc_s; var step3Model = new UnitModel(step3params); idUnitModelsList.Add(step3Model); - var distIdResult_step4 = DisturbanceIdentifier.CalculateDisturbanceVector + var distIdResult_step4 = DisturbanceCalculator.CalculateDisturbanceVector (dataSetStep3, step3Model, pidInputIdx, pidParams); idDisturbancesList.Add(distIdResult_step4); if (doConsoleDebugOut) @@ -656,7 +657,7 @@ private static Tuple GlobalSearchLinearPidGain(UnitDataSet da } // part 2: actual disturbance modeling var dataSet_altMISO = new UnitDataSet(dataSet); - var distIdResultAlt_MISO = DisturbanceIdentifier.CalculateDisturbanceVector + var distIdResultAlt_MISO = DisturbanceCalculator.CalculateDisturbanceVector (dataSet_altMISO, alternativeMISOModel, pidInputIdx, pidParams); if (false) { @@ -981,7 +982,7 @@ private static Tuple EstimateSISOdisturbanceForP int pidInputIdx_NewSystem = 0; // pidInputIdx=0 always for the new SISO system: - var candidateModelDisturbance = DisturbanceIdentifier.CalculateDisturbanceVector + var candidateModelDisturbance = DisturbanceCalculator.CalculateDisturbanceVector (dataSet_SISO, candidateModel, pidInputIdx_NewSystem, pidParams); var d_est = candidateModelDisturbance.d_est; diff --git a/Dynamic/Identification/DisturbanceIdentifier.cs b/Dynamic/Identification/DisturbanceCalculator.cs similarity index 97% rename from Dynamic/Identification/DisturbanceIdentifier.cs rename to Dynamic/Identification/DisturbanceCalculator.cs index 2cbfe78..127b091 100644 --- a/Dynamic/Identification/DisturbanceIdentifier.cs +++ b/Dynamic/Identification/DisturbanceCalculator.cs @@ -24,12 +24,12 @@ namespace TimeSeriesAnalysis.Dynamic /// /// Enum describing reasons that a disturbance may be set by logic to exactly zero /// - public enum DisturbanceSetToZeroReason + public enum DisturbanceEstimationError { /// /// default, when disturbance estimation has not run yet /// - NotRunYet=0, + None=0, /// /// Set disturbance to zero beacause simulator did not run /// @@ -58,7 +58,7 @@ public class DisturbanceIdResult /// /// The reason for setting the disturbance to zero /// - public DisturbanceSetToZeroReason zeroReason; + public DisturbanceEstimationError ErrorType; /// /// /// @@ -88,6 +88,7 @@ public class DisturbanceIdResult public DisturbanceIdResult(UnitDataSet dataSet) { N = dataSet.GetNumDataPoints(); + ErrorType = DisturbanceEstimationError.None; SetToZero(); } @@ -109,10 +110,9 @@ public void SetToZero() } /// - /// An algorithm that attempts to re-create the additive output disturbance acting on - /// a signal Y while PID-control attempts to counter-act the disturbance by adjusting its manipulated output u. + /// For a given process model and dataset, calculate the disturbance vector d by subtracting y_proc from y_meas /// - public class DisturbanceIdentifier + public class DisturbanceCalculator { /// @@ -293,7 +293,7 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat if (y_proc == null) { - result.zeroReason = DisturbanceSetToZeroReason.UnitSimulatorUnableToRun; + result.ErrorType = DisturbanceEstimationError.UnitSimulatorUnableToRun; return result; } diff --git a/Dynamic/PlantSimulator/PlantSimulator.cs b/Dynamic/PlantSimulator/PlantSimulator.cs index 527821a..fffd24f 100644 --- a/Dynamic/PlantSimulator/PlantSimulator.cs +++ b/Dynamic/PlantSimulator/PlantSimulator.cs @@ -543,6 +543,11 @@ public static (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl /// /// Simulate single model based on a unit data set + /// + /// This is a convenience function that creates a TimeSeriesDataSet, sets default names in the model and dataset that match based on unitDataset + /// The output is returned directly. + /// + /// Optionally, the result can be written to y_meas or y_sim in unitdata. /// /// contains a unit data set that must have U filled, Y_sim will be written here /// model to simulate @@ -550,11 +555,12 @@ public static (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl /// if writing to Ymeas, it is possible to add noise of the given amplitude to signal /// if true, the Y_sim of unitData has the simulation result written two i /// the seed value of the noise to be added - /// a tuple, first aa true if able to simulate, otherwise false, second is the simulated time-series - static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatableModel model,bool writeToYmeas= false, + /// a tuple, first aa true if able to simulate, otherwise false, second is the simulated time-series "y_proc" without any additive + static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatableModel model, bool writeToYmeas= false, double noiseAmplitude=0, bool addSimToUnitData=false, int seedNr=123) { + const string defaultOutputName = "output"; var inputData = new TimeSeriesDataSet(); var singleModelName = "SimulateSingle"; var modelCopy = model.Clone(singleModelName); @@ -573,18 +579,42 @@ static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl inputData.Add(uName, unitData.U.GetColumn(colIdx)); uNames.Add(uName); } - modelCopy.SetInputIDs(uNames.ToArray()); - if (modelCopy.GetOutputID() == null) - modelCopy.SetOutputID("output"); - - var sim = new PlantSimulator(new List { modelCopy }); - var isOk = sim.SimulateSingle(inputData, singleModelName, false, out var simData); + modelCopy.SetInputIDs(uNames.ToArray()); + { + inputData.Add(defaultOutputName, unitData.Y_meas); + modelCopy.SetOutputID(defaultOutputName); + } - // var isOk = PlantSimulator.SimulateSingle(inputData, modelCopy, out var simData); + var sim = new PlantSimulator(new List { modelCopy }); + var isOk = sim.SimulateSingle(inputData, singleModelName, false, out var simData); if(!isOk) return (false, null); - double[] y_sim = simData.GetValues(singleModelName, SignalType.Output_Y); + + double[] y_proc = null; + double[] y_sim = null; + + y_sim = simData.GetValues(defaultOutputName); + + if (simData.ContainsSignal(singleModelName)) + { + y_proc = simData.GetValues(singleModelName); + } + else + { + y_proc = y_sim; + } + + + // if the input included a "additive distubance" signal,there will be a internal process output + /* double[] y_proc = simData.GetValues(singleModelName); + double[] y_sim = y_proc; + + if (y_proc == null) + { + y_proc = simData.GetValues(defaultOutputName); + y_sim = simData.GetValues(singleModelName, SignalType.Output_Y); + }*/ if (noiseAmplitude > 0) { // use a specific seed here, to avoid potential issues with "random unit tests" and not-repeatable @@ -607,7 +637,7 @@ static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl unitData.Y_sim = y_sim; } } - return (isOk, y_sim); + return (isOk, y_proc); } diff --git a/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs b/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs index ed47ed1..b57eec1 100644 --- a/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs @@ -84,7 +84,8 @@ public void Static_StepDisturbance_EstimatesOk(double stepAmplitude) [TestCase(5)] public void Static_SinusDisturbance_EstimatesOk(double stepAmplitude) { - var trueDisturbance = TimeSeriesCreator.Sinus(stepAmplitude, timeBase_s*15, timeBase_s,N ); + var sinusPeriod = timeBase_s * 15; + var trueDisturbance = TimeSeriesCreator.Sinus(stepAmplitude, sinusPeriod, timeBase_s,N ); GenericDisturbanceTest(new UnitModel(staticModelParameters, "StaticProcess"), trueDisturbance); } @@ -149,7 +150,7 @@ public void GenericDisturbanceTest (UnitModel processModel, double[] trueDistur Assert.IsTrue(isOk); Assert.IsTrue(simData.ContainsSignal(processModel.GetID()),"simulated dataset should include internal process model output (pre-disturbance)"); var pidDataSet = plantSim.GetUnitDataSetForPID(inputData.Combine(simData), pidModel1); - var result = DisturbanceIdentifier.CalculateDisturbanceVector(pidDataSet, processModel); + var result = DisturbanceCalculator.CalculateDisturbanceVector(pidDataSet, processModel); if (doAssertResult) diff --git a/docs/articles/processsimulator.md b/docs/articles/processsimulator.md index f68c614..733ee79 100644 --- a/docs/articles/processsimulator.md +++ b/docs/articles/processsimulator.md @@ -55,7 +55,7 @@ Each model has a number of inputs that travel through the model, each with an id In addition, models support adding signals directly to the output. This is a feature intended for modeling disturbances. The IDs of such signals are referred to as *additive input IDs*. -### Estimating disturbances and simulating their effects +### Estimating disturbances and simulating their effects using PlantSimulator.SimulateSingle() If the ``inputData`` given to the PlantSimulator includes measured pid-outputs ``u`` and process outputs ``y`` of feedback loops, then PlantSimulator uses that