From 7b410f940bbf35e96523c41a92417c591ff318a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steinar=20Elgs=C3=A6ter?= Date: Mon, 28 Oct 2024 15:33:51 +0100 Subject: [PATCH] - added time-delay estimation to gainschedidentifier and added new unit tests to verify it - always explcitly set the seed when adding noise to a simulated Y in unit-tests to avoid race-condition issues --- Dynamic/Identification/GainSchedIdentifier.cs | 191 ++++++++++++------ .../SimulatableModels/GainSchedParameters.cs | 33 ++- .../Tests/GainSchedIdentifyTests.cs | 105 +++++++++- .../Tests/PidIdentUnitTests.cs | 8 +- TimeSeriesAnalysis/TimeSeriesDataSet.cs | 12 +- 5 files changed, 265 insertions(+), 84 deletions(-) diff --git a/Dynamic/Identification/GainSchedIdentifier.cs b/Dynamic/Identification/GainSchedIdentifier.cs index c294fc24..fb1f5b40 100644 --- a/Dynamic/Identification/GainSchedIdentifier.cs +++ b/Dynamic/Identification/GainSchedIdentifier.cs @@ -53,8 +53,8 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting UnitDataSet dataSetNoGainSched = new UnitDataSet(dataSet); GainSchedParameters GSp_noGainSchedReference = new GainSchedParameters(); GSp_noGainSchedReference.LinearGainThresholds = new double[] { }; ; - UnitModel unitModel = UnitIdentifier.IdentifyLinear(ref dataSetNoGainSched, null, false);// Todo:consider modelling with nonlinear model? - UnitParameters untiParams = unitModel.GetModelParameters(); + UnitModel unitModel = UnitIdentifier.IdentifyLinear(ref dataSetNoGainSched, null,false);// Todo:consider modelling with nonlinear model? + UnitParameters unitParams = unitModel.GetModelParameters(); bool dodebug = false; if (dodebug) { @@ -83,26 +83,28 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting } // Linear gains List GS_LinearGains1 = new List(); - GS_LinearGains1.Add(untiParams.LinearGains); + GS_LinearGains1.Add(unitParams.LinearGains); GSp_noGainSchedReference.LinearGains = GS_LinearGains1; // Time constants - double[] GS_TimeConstants_s1 = new double[] { untiParams.TimeConstant_s }; + double[] GS_TimeConstants_s1 = new double[] { unitParams.TimeConstant_s }; GSp_noGainSchedReference.TimeConstant_s = GS_TimeConstants_s1; + // time delays + GSp_noGainSchedReference.TimeDelay_s = unitParams.TimeDelay_s; + // Time constants thresholds double[] GS_TimeConstantThreshold1 = new double[] { }; GSp_noGainSchedReference.TimeConstantThresholds = GS_TimeConstantThreshold1; - GSp_noGainSchedReference.Fitting = untiParams.Fitting; + GSp_noGainSchedReference.Fitting = unitParams.Fitting; // GSp_noGainSchedReference.OperatingPoint_U = 0; - GSp_noGainSchedReference.OperatingPoint_Y = untiParams.Bias; + GSp_noGainSchedReference.OperatingPoint_Y = unitParams.Bias; ////// allYSimList.Add(dataSetNoGainSched.Y_sim); allGainSchedParams.Add(GSp_noGainSchedReference); } //////////////////////////////////////////////////// - /// // note that this is a fairly computationally heavy call // higher values means higher accuracy but at the cost of more computations.. const int globalSearchIterationsPass1 = 40;//should be above a threshold, but above that more is not better. @@ -122,20 +124,90 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting // pass 2: const bool DO_PASS2 = true; const int pass2Width = 0;//0,1 or 2, design parameter about how wide to do pass 2 aroudn pass 1 result.(higher number is at the expense of accuracy) - if (bestModelIdx_pass1 > 1 + pass2Width && bestModelIdx_pass1 < allGainSchedParams.Count() - pass2Width && DO_PASS2) + if (bestModelIdx_pass1 > 1 + pass2Width && bestModelIdx_pass1 < allGainSchedParams.Count() - pass2Width && DO_PASS2) { - double? gsSearchMin_pass2 = allGainSchedParams.ElementAt(Math.Max(bestModelIdx_pass1 - 1- pass2Width,0)).LinearGainThresholds.First(); - double? gsSearchMax_pass2 = allGainSchedParams.ElementAt(Math.Min(bestModelIdx_pass1 + 1+ pass2Width, allGainSchedParams.Count()-1)).LinearGainThresholds.First(); + double? gsSearchMin_pass2 = allGainSchedParams.ElementAt(Math.Max(bestModelIdx_pass1 - 1 - pass2Width, 0)).LinearGainThresholds.First(); + double? gsSearchMax_pass2 = allGainSchedParams.ElementAt(Math.Min(bestModelIdx_pass1 + 1 + pass2Width, allGainSchedParams.Count() - 1)).LinearGainThresholds.First(); (List potentialGainschedParametersList_pass2, List potentialYsimList_pass2) = - IdentifyGainScheduledGainsAndSingleThreshold(dataSet, gainSchedInputIndex, true, globalSearchIterationsPass2,gsSearchMin_pass2, gsSearchMax_pass2); - (var bestModel_pass2, var bestModelIdx_pass2) = ChooseBestGainScheduledModel(potentialGainschedParametersList_pass2, + IdentifyGainScheduledGainsAndSingleThreshold(dataSet, gainSchedInputIndex, true, globalSearchIterationsPass2, gsSearchMin_pass2, gsSearchMax_pass2); + (var bestModel_pass2, var bestModelIdx_pass2) = ChooseBestGainScheduledModel(potentialGainschedParametersList_pass2, potentialYsimList_pass2, ref dataSet); + EstimateTimeDelay(ref bestModel_pass2, ref dataSet); return bestModel_pass2; } else + { + EstimateTimeDelay(ref bestModel_pass1, ref dataSet); return bestModel_pass1; + } } + /// + /// Updates the given model to include a time-delay term if applicable. + /// Method is based on trial-and-error reducing the estimated time-constants while increasint the time-delay and seeing + /// if this improves the fit of the model + /// + /// estiamted paramters from other estimation, time-constants should be estiamted, but time-delays are zero, both time-delays and time-constants are updated by this method + /// + private static void EstimateTimeDelay(ref GainSchedParameters gsParams, ref UnitDataSet dataSet) + { + var minTc = (new Vec()).Min(gsParams.TimeConstant_s); + var maxTc = (new Vec()).Max(gsParams.TimeConstant_s); + var timeBase_s = dataSet.GetTimeBase(); + var vec = new Vec(); + + var resultDict = new Dictionary(); + double smallestObjFun = double.PositiveInfinity; + { + var objFun = vec.Sum(vec.Abs(vec.Subtract(dataSet.Y_sim, dataSet.Y_meas))).Value; + resultDict.Add(0, objFun); + smallestObjFun = objFun; + } + // simulate the model while increasing the time delay while decreasing the time-constants and see if this improves fit. + double timedelay_best_s = 0; + double[] y_sim_best = new double[1]; + for (double timedelay_s = timeBase_s; timedelay_s < minTc; timedelay_s += timeBase_s) + { + var copiedGsParams = new GainSchedParameters(gsParams); + copiedGsParams.TimeDelay_s = timedelay_s; + copiedGsParams.TimeConstant_s = vec.Subtract(gsParams.TimeConstant_s, timedelay_s); + + var gsIdentModel = new GainSchedModel(copiedGsParams, "ident_model"); + var identModelSim = new PlantSimulator(new List { gsIdentModel }); + var inputDataIdent = new TimeSeriesDataSet(); + + int index = 0; + foreach (var id in gsIdentModel.GetModelInputIDs()) + { + inputDataIdent.Add(identModelSim.AddExternalSignal(gsIdentModel, SignalType.External_U, index), dataSet.U.GetColumn(index)); + index++; + } + + inputDataIdent.CreateTimestamps(dataSet.GetTimeBase()); + var isOk = identModelSim.Simulate(inputDataIdent, out TimeSeriesDataSet identModelSimData); + if (isOk) + { + var currentSimY = identModelSimData.GetValues(gsIdentModel.ID, SignalType.Output_Y); + var objFun = vec.Sum(vec.Abs(vec.Subtract(currentSimY, dataSet.Y_meas))).Value; + resultDict.Add(timedelay_s, objFun); + if (objFun < smallestObjFun) + { + smallestObjFun = objFun; + timedelay_best_s = timedelay_s; + y_sim_best = currentSimY; + } + } + } + + if (timedelay_best_s > 0) + { + dataSet.Y_sim = y_sim_best; + gsParams.TimeDelay_s = timedelay_best_s; + gsParams.TimeConstant_s = vec.Subtract(gsParams.TimeConstant_s, timedelay_best_s); + } + } + + /// /// Identify a model when a given set of thresholds is already given in the supplied gsFittingSpecs /// @@ -144,13 +216,14 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting /// /// tuning data set /// the object in which the thresholds to be used are given + /// set to false to disable estimation of time delays, default is true /// - public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet, GainSchedFittingSpecs gsFittingSpecs) + public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet, GainSchedFittingSpecs gsFittingSpecs, bool doTimeDelayEstimation = true) { var vec = new Vec(); - GainSchedParameters ret = new GainSchedParameters(); - ret.GainSchedParameterIndex = gsFittingSpecs.uGainScheduledInputIndex; - ret.LinearGainThresholds = gsFittingSpecs.uGainThresholds; + GainSchedParameters retModel = new GainSchedParameters(); + retModel.GainSchedParameterIndex = gsFittingSpecs.uGainScheduledInputIndex; + retModel.LinearGainThresholds = gsFittingSpecs.uGainThresholds; // for this to work roubustly, the training set for fitting each model may need to require adding in a "span" of neighboring models. int numberOfInputs = dataSet.U.GetLength(1); double gsVarMinU, gsVarMaxU; @@ -160,7 +233,7 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet var warningNotEnoughExitationBetweenAllThresholds = false; var dataSetCopy = new UnitDataSet(dataSet); // estimate each of the gains one by one - for (int curGainIdx = 0; curGainIdx < ret.LinearGainThresholds.Count()+1; curGainIdx++) + for (int curGainIdx = 0; curGainIdx < retModel.LinearGainThresholds.Count()+1; curGainIdx++) { double[] uMinFit = new double[numberOfInputs]; double[] uMaxFit = new double[numberOfInputs]; @@ -169,23 +242,23 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet if (curGainIdx == 0) { - gsVarMinU = vec.Min(Array2D.GetColumn(dataSetCopy.U, ret.GainSchedParameterIndex)); + gsVarMinU = vec.Min(Array2D.GetColumn(dataSetCopy.U, retModel.GainSchedParameterIndex)); } else { - gsVarMinU = ret.LinearGainThresholds[curGainIdx-1]; + gsVarMinU = retModel.LinearGainThresholds[curGainIdx-1]; } - if (curGainIdx == ret.LinearGainThresholds.Count()) + if (curGainIdx == retModel.LinearGainThresholds.Count()) { - gsVarMaxU = vec.Max(Array2D.GetColumn(dataSetCopy.U, ret.GainSchedParameterIndex)); + gsVarMaxU = vec.Max(Array2D.GetColumn(dataSetCopy.U, retModel.GainSchedParameterIndex)); } else { - gsVarMaxU = ret.LinearGainThresholds[curGainIdx]; + gsVarMaxU = retModel.LinearGainThresholds[curGainIdx]; } for (int idx = 0; idx < numberOfInputs; idx++) { - if (idx == ret.GainSchedParameterIndex) + if (idx == retModel.GainSchedParameterIndex) { uMinFit[idx] = gsVarMinU; uMaxFit[idx] = gsVarMaxU; @@ -196,21 +269,16 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet uMaxFit[idx] = double.NaN; } } - var idResults = IdentifySingleLinearModelForGivenThresholds(ref dataSetCopy, uMinFit, uMaxFit, ret.GainSchedParameterIndex,u0); + var idResults = IdentifySingleLinearModelForGivenThresholds(ref dataSetCopy, uMinFit, uMaxFit, retModel.GainSchedParameterIndex,u0,false); if (idResults.NotEnoughExitationBetweenAllThresholds) warningNotEnoughExitationBetweenAllThresholds = true; - var uInsideUMinUMax = Vec.GetValuesAtIndices(dataSetCopy.U.GetColumn(ret.GainSchedParameterIndex), + var uInsideUMinUMax = Vec.GetValuesAtIndices(dataSetCopy.U.GetColumn(retModel.GainSchedParameterIndex), Index.InverseIndices(dataSetCopy.GetNumDataPoints(),dataSetCopy.IndicesToIgnore)); var uMaxObserved = (new Vec()).Max(uInsideUMinUMax); var uMinObserved = (new Vec()).Min(uInsideUMinUMax); if (idResults.LinearGains == null) { allIdsOk = false; - // Console.WriteLine("min" + gsVarMinU + " max:" + gsVarMaxU + " gain: FAILED!!");//todo:comment out when working - } - else - { - // Console.WriteLine("min" + gsVarMinU + " max:" + gsVarMaxU + " gain:" + idResults.Item1[0] + " umin:"+ uMinObserved + " umax:"+ uMaxObserved);//todo:comment out when working } linearGains.Add(idResults.LinearGains); timeConstants.Add(idResults.TimeConstant_s); @@ -218,30 +286,34 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet } // final gain:above the highest threshold - ret.LinearGains = linearGains; + retModel.LinearGains = linearGains; + // time delay if (gsFittingSpecs.uTimeConstantThresholds == null) { - ret.TimeConstant_s = new double[] { vec.Mean(timeConstants.ToArray()).Value }; + retModel.TimeConstant_s = new double[] { vec.Mean(timeConstants.ToArray()).Value }; } else if (Vec.Equal(gsFittingSpecs.uTimeConstantThresholds,gsFittingSpecs.uGainThresholds)) { - ret.TimeConstant_s = timeConstants.ToArray(); - ret.TimeConstantThresholds = gsFittingSpecs.uTimeConstantThresholds; + retModel.TimeConstant_s = timeConstants.ToArray(); + retModel.TimeConstantThresholds = gsFittingSpecs.uTimeConstantThresholds; } else { - ret.TimeConstant_s = null;//TODO: currently not supported to find timeconstants separately from gain thresholds. + retModel.TimeConstant_s = null;//TODO: currently not supported to find timeconstants separately from gain thresholds. } - ret.Fitting = new FittingInfo(); - ret.Fitting.WasAbleToIdentify = allIdsOk; + retModel.Fitting = new FittingInfo(); + retModel.Fitting.WasAbleToIdentify = allIdsOk; if (!allIdsOk) - ret.AddWarning(GainSchedIdentWarnings.UnableToIdentifySomeSubmodels); + retModel.AddWarning(GainSchedIdentWarnings.UnableToIdentifySomeSubmodels); if(warningNotEnoughExitationBetweenAllThresholds) - ret.AddWarning(GainSchedIdentWarnings.InsufficientExcitationBetweenEachThresholdToBeCertainOfGains); + retModel.AddWarning(GainSchedIdentWarnings.InsufficientExcitationBetweenEachThresholdToBeCertainOfGains); // simulate the model and determine the optimal bias term: - DetermineOperatingPointAndSimulate(ref ret, ref dataSet); - return ret; + DetermineOperatingPointAndSimulate(ref retModel, ref dataSet); + + EstimateTimeDelay(ref retModel, ref dataSet); + + return retModel; } /// @@ -257,8 +329,6 @@ private static bool DetermineOperatingPointAndSimulate(ref GainSchedParameters g var identModelSim = new PlantSimulator(new List { gsIdentModel }); var inputDataIdent = new TimeSeriesDataSet(); - //inputDataIdent.Add(identModelSim.AddExternalSignal(gsIdentModel, SignalType.External_U, (int)INDEX.FIRST), dataSet.U.GetColumn(0));// todo row 0 is not general - var vec = new Vec(); int index = 0; @@ -309,20 +379,18 @@ private static bool DetermineOperatingPointAndSimulate(ref GainSchedParameters g /// /// /// determine the operating point/bias (true by default) + /// number of global search iterations + /// minimum of the global search(used for pass 2+) + /// maximum of the global search(used for pass 2+) /// a tuple of: a) A list of all the candiate gain-scheduling parameters that have been found by global search, /// but note that these models do not include operating point, and /// b) the simulated output of each of the gain-scheduled paramters in the first list. private static (List, List) IdentifyGainScheduledGainsAndSingleThreshold(UnitDataSet dataSet, - int gainSchedInputIndex, bool estimateOperatingPoint=true, int globalSearchIterations = 40, double? gsSearchMin=null, double? gsSearchMax= null ) + int gainSchedInputIndex, bool estimateOperatingPoint=true, + int globalSearchIterations = 40, double? gsSearchMin=null, double? gsSearchMax= null ) { - - // TODO: since there are magic number below, this method may in some cases not work as intended, ideally there - // should be no magic numbers. - - // the number of thresholds to try between the minium and maximum - // const int globalSearchIterations = 40; // how big a span of the range of u in the dataset to span with trehsolds (should be <= 1.00 and >0.) - // usually it is pointless to search for thresholds at the edge of the dataset, so should be <1 + // usually it is pointless to search for thresholds at the edge of the dataset, so should be <1, but always much larger than 0. const double GS_RANGE_SEARCH_FRAC = 0.70; // when didving the dataset into "a" and "b" above and below the threshold, how many percentage overlap @@ -438,9 +506,10 @@ private static (List, List) IdentifyGainScheduled /// the upper threshold of the gain-scheduled input /// index of the variable in the model that is gain-scheduled. by the fault 0 /// if set to null, the u between uLowerThreshold and uHigherThreshold is used + /// if set to true, identification also considers time delay (increases computational load, use sparingly) /// return a tuple, the first linear gains is null if unable to identify private static GainSchedSubModelResults IdentifySingleLinearModelForGivenThresholds(ref UnitDataSet dataSet,double[] uLowerThreshold, double[] uHigherThreshold, - int gainSchedVarIndex, double? u0= null ) + int gainSchedVarIndex, double? u0= null, bool doTimeDelayEstimation = false ) { var results = new GainSchedSubModelResults(); @@ -456,12 +525,21 @@ private static GainSchedSubModelResults IdentifySingleLinearModelForGivenThresho results.NotEnoughExitationBetweenAllThresholds = false; + int[] incomingIndToIgnore = null; + if (dataSet.IndicesToIgnore != null) + { + incomingIndToIgnore = new int[dataSet.IndicesToIgnore.Count]; + dataSet.IndicesToIgnore.CopyTo( incomingIndToIgnore); + } // if there is close to as many thresholds as there are steps in data, then there may not actually be enough data in the dataset to // this while loop looks at the range betwen uMin and uMax for the data found between the upper and lower thresold. If the range // is low, then the range of u to fit against is attempted increased to capture information as u steps between thresholds. while (uRangeInChosenDataset < MIN_uRangeInChosenDataset && tuningSetSpanOutsideThreshold_prc < MAX_tuningSetSpanOutsideThreshold_prc) { - dataSet.IndicesToIgnore = null;// TODO: store any incoming? + if (incomingIndToIgnore == null) + dataSet.IndicesToIgnore = null; + else + dataSet.IndicesToIgnore = new List(incomingIndToIgnore); fittingSpecs.U_min_fit = uLowerThreshold; fittingSpecs.U_max_fit = uHigherThreshold; fittingSpecs.U_min_fit[gainSchedVarIndex] = uLowerThreshold[gainSchedVarIndex] - fitRange * tuningSetSpanOutsideThreshold_prc/100; @@ -481,15 +559,12 @@ private static GainSchedSubModelResults IdentifySingleLinearModelForGivenThresho if (tuningSetSpanOutsideThreshold_prc > 0) { results.NotEnoughExitationBetweenAllThresholds = true; - // Console.WriteLine("uRangeInChosenDataset:" + uRangeInChosenDataset + " tuningSetSpanOutsideThreshold_prc:" + tuningSetSpanOutsideThreshold_prc); } tuningSetSpanOutsideThreshold_prc += STEP_tuningSetSpanOutsideThreshold_prc; } - // NB! save a lot of computational time by not doing time delay estimation here! - // var unitModel = UnitIdentifier.Identify(ref dataSet, fittingSpecs, false); - bool doTimeDelayEstimation = false; - var unitModel = UnitIdentifier.IdentifyLinear(ref dataSet, fittingSpecs, doTimeDelayEstimation); + // NB! the algorithm saves a lot of computational time by not doing time delay estimation here! + var unitModel = UnitIdentifier.IdentifyLinear(ref dataSet, fittingSpecs, false); bool doDebug = false; // should be false unless debugging if (doDebug) diff --git a/Dynamic/SimulatableModels/GainSchedParameters.cs b/Dynamic/SimulatableModels/GainSchedParameters.cs index e7eb61d7..3bb894cd 100644 --- a/Dynamic/SimulatableModels/GainSchedParameters.cs +++ b/Dynamic/SimulatableModels/GainSchedParameters.cs @@ -68,17 +68,6 @@ public class GainSchedParameters : ModelParametersBaseClass /// public List LinearGainUnc { get; set; } = null; - /// - /// The working point of the model, the value of each U around which the model is localized. - /// If value is nullc> then no U0 is used in the model - /// - // public double[] U0 { get; set; } = null; // TODO: consider removing for gain-scheduled model... - - - /// - /// Note that for gain-scheduled models this is a private paramter that is calcualted - /// - // private double Bias = 0; /// /// The "operating point" specifies the value y that the model should have for the gain scheduled input u. @@ -96,6 +85,28 @@ public GainSchedParameters() errorsAndWarningMessages = new List(); } + public GainSchedParameters(GainSchedParameters existingModel) + { + Y_min = existingModel.Y_min; + Y_max = existingModel.Y_max; + TimeConstantThresholds = existingModel.TimeConstantThresholds; + TimeConstant_s = existingModel.TimeConstant_s; + TimeConstantUnc_s = existingModel.TimeConstantUnc_s; + LinearGains = existingModel.LinearGains; + LinearGainUnc = existingModel.LinearGainUnc; + LinearGainThresholds = existingModel.LinearGainThresholds; + Fitting = existingModel.Fitting; + errorsAndWarningMessages = existingModel.errorsAndWarningMessages; + OperatingPoint_U = existingModel.OperatingPoint_U; + OperatingPoint_Y = existingModel.OperatingPoint_Y; + GainSchedParameterIndex = existingModel.GainSchedParameterIndex; + TimeDelay_s = existingModel.TimeDelay_s; + FittingSpecs = existingModel.FittingSpecs; + + + } + + /// /// Returns the bias calculated from OperatingPoint_U, OperatingPoint_Y; /// diff --git a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs b/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs index ba85271f..d6779273 100644 --- a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs @@ -11,6 +11,8 @@ using System.Reflection; using Accord.Math.Transforms; +using System.Globalization; + namespace TimeSeriesAnalysis.Test.SysID { public class GainSchedIdentifyTests @@ -18,7 +20,29 @@ public class GainSchedIdentifyTests const int timeBase_s = 1; const double TimeConstantAllowedDev_s = 3.5; + private void ConoleOutResult(GainSchedParameters trueParams, GainSchedParameters estParams) + { + string delimTxt = " vs true value: "; + string accuracy = "F2"; + + for (int i = 0; i < estParams.TimeConstant_s.Length; i++) + { + Console.Write("Est timeconstant "+(i+1)+":" + estParams.TimeConstant_s.ElementAt(i).ToString(accuracy, CultureInfo.InvariantCulture) + + delimTxt + trueParams.TimeConstant_s.ElementAt(i).ToString(accuracy, CultureInfo.InvariantCulture) + "\r\n"); + } + for (int i = 0; i < estParams.LinearGains.Count; i++) + { + Console.Write("Est gain " + (i + 1) + ":" + estParams.LinearGains.ElementAt(i).ElementAt(0).ToString(accuracy, CultureInfo.InvariantCulture) + + delimTxt + estParams.LinearGains.ElementAt(i).ElementAt(0).ToString(accuracy, CultureInfo.InvariantCulture) + "\r\n"); + } + Console.Write("Op point Y:" + estParams.OperatingPoint_Y.ToString(accuracy, CultureInfo.InvariantCulture) + + delimTxt + trueParams.OperatingPoint_Y.ToString(accuracy, CultureInfo.InvariantCulture) + "\r\n"); + Console.Write("threshold :" + estParams.LinearGainThresholds.First().ToString(accuracy) + + delimTxt + trueParams.LinearGainThresholds.First().ToString(accuracy, CultureInfo.InvariantCulture) + "\r\n"); + Console.Write("time-delay :" + estParams.TimeDelay_s.ToString(accuracy, CultureInfo.InvariantCulture) + + delimTxt + trueParams.TimeDelay_s.ToString(accuracy, CultureInfo.InvariantCulture) + "\r\n"); + } // note that the tolerance seems to be linear with the noise in the data // five varying gains @@ -43,7 +67,7 @@ public void FiveGains_CorrectGainsReturned(int ver, int expectedNumWarnings, dou LinearGains = new List { new double[] { 0.5 }, new double[] { 1 }, new double[] { 3 }, new double[] { 4.5 }, new double[] { 6 }, new double[] { 9 } }, LinearGainThresholds = new double[] { 2.5, 4.5, 6.5, 8.5, 10.5 }, TimeDelay_s = 0, - OperatingPoint_U = 5, + OperatingPoint_U = 0, OperatingPoint_Y = 4, GainSchedParameterIndex = 0 }; @@ -57,7 +81,7 @@ public void FiveGains_CorrectGainsReturned(int ver, int expectedNumWarnings, dou LinearGains = new List { new double[] { 2 }, new double[] { 2 }, new double[] { 2}, new double[] { 2 }, new double[] { 2 }, new double[] { 2 } }, LinearGainThresholds = new double[] { 2.5, 4.5, 6.5, 8.5, 10.5 }, TimeDelay_s = 0, - OperatingPoint_U = 5, + OperatingPoint_U = 0, OperatingPoint_Y = 4, GainSchedParameterIndex = 0 }; @@ -78,7 +102,7 @@ public void FiveGains_CorrectGainsReturned(int ver, int expectedNumWarnings, dou // Act var isSimulatable = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); - simData.AddNoiseToSignal(SignalNamer.GetSignalName(refModel.ID, SignalType.Output_Y, 0),noiseAmplitude); + simData.AddNoiseToSignal(SignalNamer.GetSignalName(refModel.ID, SignalType.Output_Y, 0),noiseAmplitude,123); Assert.IsTrue(isSimulatable); var dataSet = new UnitDataSet(); @@ -133,6 +157,79 @@ private void DiffLessThan(double trueVal, double testVal, double tolerancePrc,in } } + [TestCase(3)] + [TestCase(5)] + [TestCase(7)] + + public void TimeDelay_TDEstOk(int timeDelaySamples) + { + double noiseAmp = 0.25; + int N = 300; + // Arrange + var unitData = new UnitDataSet("test"); + double[] u1 = TimeSeriesCreator.ThreeSteps(N / 5, N / 3, N / 2, N, -2, -1, 0, 1); + double[] u2 = TimeSeriesCreator.ThreeSteps(3 * N / 5, 2 * N / 3, 4 * N / 5, N, 0, 1, 2, 3); + double[] u = u1.Zip(u2, (x, y) => x + y).ToArray(); + double[,] U = Array2D.CreateFromList(new List { u }); + unitData.U = U; + unitData.Times = TimeSeriesCreator.CreateDateStampArray( + new DateTime(2000, 1, 1), timeBase_s, N); + + double threshold =2; + + //reference model + GainSchedParameters trueGSparams = new GainSchedParameters + { + TimeConstant_s = new double[] { 3, 10 }, + TimeConstantThresholds = new double[] { threshold }, + LinearGains = new List { new double[] { -2 }, new double[] { 3 } }, + LinearGainThresholds = new double[] { threshold }, + TimeDelay_s = timeBase_s* timeDelaySamples, + }; + trueGSparams.OperatingPoint_Y = -1.34; + + GainSchedModel trueModel = new GainSchedModel(trueGSparams, "Correct gain sched model"); + var correct_plantSim = new PlantSimulator(new List { trueModel }); + var inputData = new TimeSeriesDataSet(); + inputData.Add(correct_plantSim.AddExternalSignal(trueModel, SignalType.External_U, (int)INDEX.FIRST), u); + inputData.CreateTimestamps(timeBase_s); + var isOk = correct_plantSim.Simulate(inputData, out TimeSeriesDataSet refSimData); + SISOTests.CommonAsserts(inputData, refSimData, correct_plantSim); + double[] simY1 = refSimData.GetValues(trueModel.GetID(), SignalType.Output_Y); + unitData.Y_meas = (new Vec()).Add(Vec.Rand(simY1.Length, -noiseAmp, noiseAmp, 454), simY1); + + // Act + GainSchedParameters idParams = GainSchedIdentifier.Identify(unitData); + + // plot + bool doPlot = false; + if (doPlot) + { + Shared.EnablePlots(); + Plot.FromList(new List { + unitData.Y_meas , + unitData.Y_sim, + unitData.U.GetColumn(0) }, + new List { "y1=y_meas", "y1=y_ident", "y3=u1" }, + timeBase_s, + "GainSched - timeconstant - "); + Shared.DisablePlots(); + } + + ConoleOutResult(trueGSparams, idParams); + + // Assert + /* int min_number_of_gains = Math.Min(idParams.LinearGainThresholds.Length, trueGSparams.LinearGainThresholds.Length); + for (int k = 0; k < min_number_of_gains; k++) + { + Console.WriteLine("identified threshold: " + idParams.LinearGainThresholds[k].ToString("F3") + "true threshold: " + trueGSparams.LinearGainThresholds[k].ToString("F3")); + Assert.That(Math.Abs(idParams.LinearGainThresholds[k] - trueGSparams.LinearGainThresholds[k]), Is.LessThanOrEqualTo(linearGainTresholdTol), + "There are too large differences in the linear gain threshold " + k.ToString()); + }*/ + } + + + // note that the input varies from -2 to 4 here, so threshold beyond that are not identifiable, and at the edges they are also hard to identify. [TestCase(-0.5, 0.055)] @@ -167,6 +264,8 @@ public void TwoGainsVaryingTreshold_ThresholdEstOk(double gain_sched_threshold, LinearGainThresholds = new double[] { gain_sched_threshold }, TimeDelay_s = 0, }; + trueGSparams.OperatingPoint_Y = -1.34; + GainSchedModel trueModel = new GainSchedModel(trueGSparams, "Correct gain sched model"); var correct_plantSim = new PlantSimulator(new List { trueModel }); var inputData = new TimeSeriesDataSet(); diff --git a/TimeSeriesAnalysis.Tests/Tests/PidIdentUnitTests.cs b/TimeSeriesAnalysis.Tests/Tests/PidIdentUnitTests.cs index 8340c3ad..ee2f8fc6 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PidIdentUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/PidIdentUnitTests.cs @@ -55,7 +55,7 @@ public void SetpointStep_WNoise_KpAndTiEstimatedOk(double ySetAmplitude, double inputData.Add(processSim.AddExternalSignal(pidModel1, SignalType.Setpoint_Yset), TimeSeriesCreator.Step(N/7, N,50,50+ySetAmplitude)); inputData.CreateTimestamps(timeBase_s,t0); var isOk = processSim.Simulate(inputData,out TimeSeriesDataSet simData); - simData.AddNoiseToSignal("SubProcess1-Output_Y", yNoiseAmplitude); + simData.AddNoiseToSignal("SubProcess1-Output_Y", yNoiseAmplitude,0); Assert.IsTrue(isOk); var pidDataSet = processSim.GetUnitDataSetForPID(inputData.Combine(simData), pidModel1); @@ -107,7 +107,7 @@ public void SetpointStep_WNoise_Downsampled_KpAndTiEstimatedOk(int downsampleFac inputData.Add(processSim.AddExternalSignal(pidModel1, SignalType.Setpoint_Yset), TimeSeriesCreator.Step(N / 7, N, 50, 50 + ySetAmplitude)); inputData.CreateTimestamps(timeBase_s, t0); var isOk = processSim.Simulate(inputData, out TimeSeriesDataSet simData); - simData.AddNoiseToSignal("SubProcess1-Output_Y", yNoiseAmplitude); + simData.AddNoiseToSignal("SubProcess1-Output_Y", yNoiseAmplitude,0); Assert.IsTrue(isOk); var combinedData = inputData.Combine(simData); @@ -162,7 +162,7 @@ public void DistStep_WNoise_KpAndTiEstimatedOk(double stepAmplitude, double yNoi inputData.Add(processSim.AddExternalSignal(processModel1, SignalType.Disturbance_D), TimeSeriesCreator.Step(N/2,N,0,stepAmplitude)); inputData.CreateTimestamps(timeBase_s); var isOk = processSim.Simulate(inputData, out TimeSeriesDataSet simData); - simData.AddNoiseToSignal("SubProcess1-Output_Y", yNoiseAmplitude); + simData.AddNoiseToSignal("SubProcess1-Output_Y", yNoiseAmplitude,890978); Assert.IsTrue(isOk); var pidDataSet = processSim.GetUnitDataSetForPID(inputData.Combine(simData), pidModel1); @@ -261,7 +261,7 @@ public void DistStep_WNoise_Downsampled_KpAndTiEstimatedOk(int downsampleFactor, inputData.Add(processSim.AddExternalSignal(processModel1, SignalType.Disturbance_D), TimeSeriesCreator.Step(N / 2, N, 0, stepAmplitude)); inputData.CreateTimestamps(timeBase_s); var isOk = processSim.Simulate(inputData, out TimeSeriesDataSet simData); - simData.AddNoiseToSignal("SubProcess1-Output_Y", noiseAmplitude); + simData.AddNoiseToSignal("SubProcess1-Output_Y", noiseAmplitude,495495); var combinedData = inputData.Combine(simData); var downsampleData = combinedData.CreateDownsampledCopy(downsampleFactor); // ----do not use inputData or simData below this line---- diff --git a/TimeSeriesAnalysis/TimeSeriesDataSet.cs b/TimeSeriesAnalysis/TimeSeriesDataSet.cs index 0e4d8c77..250483e0 100644 --- a/TimeSeriesAnalysis/TimeSeriesDataSet.cs +++ b/TimeSeriesAnalysis/TimeSeriesDataSet.cs @@ -171,19 +171,15 @@ public bool AddDataPoint(string signalID, int idx, double value) /// Adds noise to a given signal in the datset. /// (This is mainly intended for testing identification algorithms against simulated data.) /// - /// - /// - /// + /// name of signal to have noise added to it + /// the amplutide of noise, the noise will be [-noiseAmplitude, noiseAmplitude] + /// a integer seed number is /// - public bool AddNoiseToSignal(string signalName, double noiseAmplitude, int? seed = null) + public bool AddNoiseToSignal(string signalName, double noiseAmplitude, int seed) { if (!dataset.ContainsKey(signalName)) return false; - if (seed == null) - { - seed = 0; - } dataset[signalName] = new Vec().Add(dataset[signalName], Vec.Rand(N.Value, -noiseAmplitude, noiseAmplitude, seed)); return true; }