From 3a2c89487a0d76212f4dde49a2fe152748097968 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steinar=20Elgs=C3=A6ter?= Date: Mon, 28 Oct 2024 11:24:01 +0100 Subject: [PATCH] -refactoring uni tests for gainschedidentifier - fine-tuning some of the pass 2 parameters in gainschedidentifier --- Dynamic/Identification/GainSchedIdentifier.cs | 41 +-- Dynamic/SimulatableModels/GainSchedModel.cs | 10 +- .../Tests/GainSchedIdentifyTests.cs | 314 +++--------------- 3 files changed, 68 insertions(+), 297 deletions(-) diff --git a/Dynamic/Identification/GainSchedIdentifier.cs b/Dynamic/Identification/GainSchedIdentifier.cs index 157d2e8..c294fc2 100644 --- a/Dynamic/Identification/GainSchedIdentifier.cs +++ b/Dynamic/Identification/GainSchedIdentifier.cs @@ -102,9 +102,14 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting 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. + const int globalSearchIterationsPass2 = 20;//should be above a threshold, but above that more is not better. + (List potentialGainschedParametersList, List potentialYsimList) = - IdentifyGainScheduledGainsAndSingleThreshold(dataSet, gainSchedInputIndex); + IdentifyGainScheduledGainsAndSingleThreshold(dataSet, gainSchedInputIndex,true, globalSearchIterationsPass1); // TODO: extend the method to consider multiple thresholds by calling the above method recursively? //////////////////////////////////////////////////// @@ -112,19 +117,19 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting allGainSchedParams.AddRange(potentialGainschedParametersList); allYSimList.AddRange(potentialYsimList); - (var bestModel_pass1, var bestModelIdx_pass1) = ChooseBestGainScheduledModel(allGainSchedParams, allYSimList, dataSet); + (var bestModel_pass1, var bestModelIdx_pass1) = ChooseBestGainScheduledModel(allGainSchedParams, allYSimList, ref dataSet); // pass 2: const bool DO_PASS2 = true; - const int pass2Width = 1;//0,1 or 2, design parameter about how wide to do pass 2 aroudn pass 1 result. + 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) { - double? gsSearchMin_pass2 = allGainSchedParams.ElementAt(bestModelIdx_pass1-1- pass2Width).LinearGainThresholds.First(); - double? gsSearchMax_pass2 = allGainSchedParams.ElementAt(bestModelIdx_pass1+1+ pass2Width).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, gsSearchMin_pass2, gsSearchMax_pass2); + IdentifyGainScheduledGainsAndSingleThreshold(dataSet, gainSchedInputIndex, true, globalSearchIterationsPass2,gsSearchMin_pass2, gsSearchMax_pass2); (var bestModel_pass2, var bestModelIdx_pass2) = ChooseBestGainScheduledModel(potentialGainschedParametersList_pass2, - potentialYsimList_pass2, dataSet); + potentialYsimList_pass2, ref dataSet); return bestModel_pass2; } else @@ -308,14 +313,14 @@ private static bool DetermineOperatingPointAndSimulate(ref GainSchedParameters g /// 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, 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 THRESHOLD_GLOBALSEARCH_MAX_IT = 40; + // 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 const double GS_RANGE_SEARCH_FRAC = 0.70; @@ -335,24 +340,24 @@ private static (List, List) IdentifyGainScheduled List potentialYSimList = new List(); // begin setting up global search - double[] potential_gainthresholds = new double[THRESHOLD_GLOBALSEARCH_MAX_IT]; + double[] potential_gainthresholds = new double[globalSearchIterations]; // default global search based on the range of values in the given dataset if (!gsSearchMin.HasValue || !gsSearchMax.HasValue) { double gsRange = gsVarMaxU - gsVarMinU; double gsSearchRange = gsRange * GS_RANGE_SEARCH_FRAC; - for (int k = 0; k < THRESHOLD_GLOBALSEARCH_MAX_IT; k++) + for (int k = 0; k < globalSearchIterations; k++) { - potential_gainthresholds[k] = gsVarMinU + (1-GS_RANGE_SEARCH_FRAC) * gsRange / 2 + gsSearchRange * ((double)k / THRESHOLD_GLOBALSEARCH_MAX_IT); + potential_gainthresholds[k] = gsVarMinU + (1-GS_RANGE_SEARCH_FRAC) * gsRange / 2 + gsSearchRange * ((double)k / globalSearchIterations); } } else // search based on the given min and max (used for second or third passes to improve accuracy) { double gsRange = gsSearchMax.Value - gsSearchMin.Value; double gsSearchRange = gsRange ; - for (int k = 0; k < THRESHOLD_GLOBALSEARCH_MAX_IT; k++) + for (int k = 0; k < globalSearchIterations; k++) { - potential_gainthresholds[k] = gsSearchMin.Value + gsSearchRange * ((double)k / THRESHOLD_GLOBALSEARCH_MAX_IT); + potential_gainthresholds[k] = gsSearchMin.Value + gsSearchRange * ((double)k / globalSearchIterations); } } @@ -523,10 +528,10 @@ private static GainSchedSubModelResults IdentifySingleLinearModelForGivenThresho /// /// all candidate gain-scheduled parameters /// list of all the simulated y corresponding to each parameter set in the first input - /// + /// the dataset to be returned, where y_sim is to be set based on the best model. /// a tuple with the paramters of the "best" model and the index of this model among the list provided static private (GainSchedParameters, int) ChooseBestGainScheduledModel(List allGainSchedParams, - List allYsimList, UnitDataSet dataSet) + List allYsimList, ref UnitDataSet dataSet) { var vec = new Vec(); GainSchedParameters BestGainSchedParams = null; @@ -539,7 +544,6 @@ static private (GainSchedParameters, int) ChooseBestGainScheduledModel(List x.ToString("F2", CultureInfo.InvariantCulture)))); Shared.DisablePlots(); - // todo: add pause? Thread.Sleep(1000); } if (rms < lowest_rms) @@ -576,8 +579,8 @@ static private (GainSchedParameters, int) ChooseBestGainScheduledModel(List= gainSchedStartModelIdx) { - // rememebr if there are N gains, there are N-1 gain tresholds + // remember if there are N gains, there are N-1 gain tresholds deltaU = (modelParameters.LinearGainThresholds[gainSchedStartModelIdx-1] - uGainSched_Start); gainsToReturn += deltaU * modelParameters.LinearGains.ElementAt(gainSchedStartModelIdx)[inputIndex]; } - // else - // { - // Console.WriteLine("wtf"); - // } // middle entire portions for (int curGainSchedModIdx = gainSchedStartModelIdx - 1; curGainSchedModIdx > gainSchedEndModelIdx; curGainSchedModIdx--) { - deltaU = (modelParameters.LinearGainThresholds[curGainSchedModIdx-1] - - modelParameters.LinearGainThresholds[curGainSchedModIdx - 2]); + deltaU = (modelParameters.LinearGainThresholds[curGainSchedModIdx] - + modelParameters.LinearGainThresholds[curGainSchedModIdx - 1]); gainsToReturn += deltaU * modelParameters.LinearGains.ElementAt(curGainSchedModIdx)[inputIndex]; } // last portions (might be a treshold to inbetween two tresholds) diff --git a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs b/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs index a691059..ba85271 100644 --- a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs @@ -9,6 +9,7 @@ using System.Diagnostics; using System; using System.Reflection; +using Accord.Math.Transforms; namespace TimeSeriesAnalysis.Test.SysID { @@ -28,7 +29,7 @@ public class GainSchedIdentifyTests [TestCase(2, 0, 1, 0.0, Description = "Two steps for every threshold(five thresholds)")] [TestCase(2, 0, 5, 1.0, Description = "Two steps for every threshold(five thresholds)")]//note here that the tolerance can be set much lower! [TestCase(2, 0, 10, 2.0, Description = "Two steps for every threshold(five thresholds)")]//note here that the tolerance can be set much lower! - public void GainEst_FiveGains_CorrectGainsReturned(int ver, int expectedNumWarnings, double gainTolerancePrc, double noiseAmplitude ) + public void FiveGains_CorrectGainsReturned(int ver, int expectedNumWarnings, double gainTolerancePrc, double noiseAmplitude ) { const int N = 100;//Note, that the actual dataset is four times this value. GainSchedParameters refParams = new GainSchedParameters(); @@ -132,90 +133,20 @@ private void DiffLessThan(double trueVal, double testVal, double tolerancePrc,in } } - [TestCase()] - public void GainAndThreshold_GainsNotLargerThanTheBiggestPossibleGain() - { - int N = 500; - - // Arrange - var unitData = new UnitDataSet("test"); /* Create an instance of TimeSeries with test data */ - double[] u1 = TimeSeriesCreator.ThreeSteps(N / 5, N / 3, N / 2, N, 0, 1, 2, 3); - 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); - - var correct_gain_sched_parameters = new GainSchedParameters - { - TimeConstant_s = new double[] { 10 }, - TimeConstantThresholds = new double[] { }, - LinearGains = new List { new double[] { 1 }, new double[] { 6 } }, - LinearGainThresholds = new double[] { 3.1 }, - TimeDelay_s = 0, - }; - GainSchedModel correct_model = new GainSchedModel(correct_gain_sched_parameters, "Correct gain sched model"); - var correct_plantSim = new PlantSimulator(new List { correct_model }); - var inputData = new TimeSeriesDataSet(); - inputData.Add(correct_plantSim.AddExternalSignal(correct_model, SignalType.External_U, (int)INDEX.FIRST), u); - inputData.CreateTimestamps(timeBase_s); - var CorrectisSimulatable = correct_plantSim.Simulate(inputData, out TimeSeriesDataSet CorrectsimData); - SISOTests.CommonAsserts(inputData, CorrectsimData, correct_plantSim); - double[] simY1 = CorrectsimData.GetValues(correct_model.GetID(), SignalType.Output_Y); - unitData.Y_meas = simY1; - - // Act - var best_params = GainSchedIdentifier.Identify(unitData); - double current_abs_value = 0; - double largest_gain_amplitude = 0; - for (int k = 0; k < best_params.LinearGains.Count; k++) - { - current_abs_value = Math.Sqrt(best_params.LinearGains[k][0] * best_params.LinearGains[k][0]); - if (current_abs_value > largest_gain_amplitude) - { - largest_gain_amplitude = current_abs_value; - } - } - - double largest_correct_gain_amplitude = 0; - for (int k = 0; k < correct_gain_sched_parameters.LinearGains.Count; k++) - { - current_abs_value = Math.Sqrt(correct_gain_sched_parameters.LinearGains[k][0] * correct_gain_sched_parameters.LinearGains[k][0]); - if (current_abs_value > largest_correct_gain_amplitude) - { - largest_correct_gain_amplitude = current_abs_value; - } - } + // 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)] + [TestCase(-0.2, 0.055)] + [TestCase(0.2, 0.045)] + [TestCase(0.5, 0.04)] + [TestCase(1.0, 0.01)] + [TestCase(2.0, 0.015)] + [TestCase(2.5, 0.015)] + [TestCase(3.0, 0.015) ] - bool doPlot = true; - if (doPlot) - { - Shared.EnablePlots(); - Plot.FromList(new List { - simY1, - unitData.U.GetColumn(0) }, - new List { "y1=correct_model", "y3=u1" }, - timeBase_s, - "GainSched - Max Threshold"); - Shared.DisablePlots(); - } - // Assert - Assert.That(largest_gain_amplitude, Is.LessThanOrEqualTo(largest_correct_gain_amplitude), - "The largest gain in the best fitting model cannot exceed the largest gain amplitude of the correct model"); - - } - - /* [TestCase(1, -1.5)] - [TestCase(2, -1.0)] - [TestCase(3, -0.5)] - [TestCase(4, 1.0)] - [TestCase(5, 2.5)] - [TestCase(6, 3.0)]*/ - [TestCase(7, 4.0)] - public void GainAndThreshold_LinearGainThresholdAtReasonablePlace(int ver, double gain_sched_threshold) + public void TwoGainsVaryingTreshold_ThresholdEstOk(double gain_sched_threshold, double linearGainTresholdTol ) { + double noiseAmp = 0.25; int N = 300; // Arrange var unitData = new UnitDataSet("test"); @@ -228,7 +159,7 @@ public void GainAndThreshold_LinearGainThresholdAtReasonablePlace(int ver, doubl new DateTime(2000, 1, 1), timeBase_s, N); //reference model - GainSchedParameters correct_gain_sched_parameters = new GainSchedParameters + GainSchedParameters trueGSparams = new GainSchedParameters { TimeConstant_s = new double[] { 3, 10 }, TimeConstantThresholds = new double[] { gain_sched_threshold }, @@ -236,43 +167,42 @@ public void GainAndThreshold_LinearGainThresholdAtReasonablePlace(int ver, doubl LinearGainThresholds = new double[] { gain_sched_threshold }, TimeDelay_s = 0, }; - GainSchedModel correct_model = new GainSchedModel(correct_gain_sched_parameters, "Correct gain sched model"); - var correct_plantSim = new PlantSimulator(new List { correct_model }); + 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(correct_model, SignalType.External_U, (int)INDEX.FIRST), u); + inputData.Add(correct_plantSim.AddExternalSignal(trueModel, SignalType.External_U, (int)INDEX.FIRST), u); inputData.CreateTimestamps(timeBase_s); - var CorrectisSimulatable = correct_plantSim.Simulate(inputData, out TimeSeriesDataSet CorrectsimData); - SISOTests.CommonAsserts(inputData, CorrectsimData, correct_plantSim); - double[] simY1 = CorrectsimData.GetValues(correct_model.GetID(), SignalType.Output_Y); - unitData.Y_meas = simY1; + 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, (int)Math.Ceiling(2*gain_sched_threshold+45)), simY1); // Act - GainSchedParameters best_params = GainSchedIdentifier.Identify(unitData); - GainSchedModel best_model = new GainSchedModel(best_params, "Best fitting model"); - var best_plantSim = new PlantSimulator(new List { best_model }); - inputData.Add(best_plantSim.AddExternalSignal(best_model, SignalType.External_U, (int)INDEX.FIRST), u); - var IdentifiedisSimulatable = best_plantSim.Simulate(inputData, out TimeSeriesDataSet IdentifiedsimData); - SISOTests.CommonAsserts(inputData, IdentifiedsimData, best_plantSim); - double[] simY2 = IdentifiedsimData.GetValues(best_model.GetID(), SignalType.Output_Y); - // Number of inputs can be determined from the TimeSeries object, assuming it provides a way to determine this - int numberOfInputs = unitData.U.GetNColumns(); // Example property, replace with actual implementation + 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 - Threshold - " + gain_sched_threshold.ToString("F2")); + Shared.DisablePlots(); + } + // Assert - int min_number_of_gains = Math.Min(best_params.LinearGainThresholds.Length, correct_gain_sched_parameters.LinearGainThresholds.Length); + int min_number_of_gains = Math.Min(idParams.LinearGainThresholds.Length, trueGSparams.LinearGainThresholds.Length); for (int k = 0; k < min_number_of_gains; k++) { - Assert.That(Math.Pow(best_params.LinearGainThresholds[k] - correct_gain_sched_parameters.LinearGainThresholds[k], 2), Is.LessThanOrEqualTo(0.5), + 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()); } - /* - Shared.EnablePlots(); - Plot.FromList(new List { - simY1, - simY2, - unitData.U.GetColumn(0) }, - new List { "y1=correct_model", "y1=best_model", "y3=u1" }, - timeBase_s, - "GainSched - Threshold at reasonable place - " + ver.ToString()); - Shared.DisablePlots();*/ } [TestCase(3, 10, 0, Description= "identify gain and time constants, zero bias, thresholds are GIVEN")] @@ -280,7 +210,7 @@ public void GainAndThreshold_LinearGainThresholdAtReasonablePlace(int ver, doubl [TestCase(3, 10, 2, Description = "identify gain and time constants AND THRESHOLDS, zero bias, ")] [TestCase(3, 10, 3, Description = "same as ver2, except non-zero bias")] - public void TimeConstant_TwoGains_TCAndThresholdFoundOk(double TimeConstant1_s, + public void TwoGains_TCAndThresholdFoundOk(double TimeConstant1_s, double TimeConstant2_s, int ver) { int N = 300; @@ -363,10 +293,7 @@ public void TimeConstant_TwoGains_TCAndThresholdFoundOk(double TimeConstant1_s, Console.Write("threshold :" + est_params.LinearGainThresholds.First().ToString("F2") + " vs true value: " + trueParams.LinearGainThresholds.First() + "\r\n"); - - // Asserts - Assert.IsTrue(Math.Abs(est_params.LinearGainThresholds.First() - trueParams.LinearGainThresholds.First()) < threshold_tol); SISOTests.CommonAsserts(inputData, IdentifiedsimData, estPlantSim); @@ -379,163 +306,8 @@ public void TimeConstant_TwoGains_TCAndThresholdFoundOk(double TimeConstant1_s, "Too high time constant " + k.ToString()); } Assert.That(Math.Abs(est_params.OperatingPoint_Y - trueParams.OperatingPoint_Y) < operatingy_tol); - - - } - - [TestCase(1, 1.5)] - /*[TestCase(2, 2.0)] - [TestCase(3, 2.5)] - [TestCase(4, 3.0)] - [TestCase(5, 3.5)] - [TestCase(6, 4.0)]*/ - // [TestCase(7, 4.5)] - public void GainAndThreshold_ThresholdsWithinUminAndUmax(int ver, double gain_sched_threshold) - { - int N = 250; - // Arrange - var unitData = new UnitDataSet("test"); /* Create an instance of TimeSeries with test data */ - 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); - - // reference model - GainSchedParameters correct_gain_sched_parameters = new GainSchedParameters - { - TimeConstant_s = new double[] { 3, 10 }, - TimeConstantThresholds = new double[] { gain_sched_threshold }, - LinearGains = new List { new double[] { 1 }, new double[] { 3 } }, - LinearGainThresholds = new double[] { gain_sched_threshold }, - TimeDelay_s = 0, - }; - GainSchedModel correct_model = new GainSchedModel(correct_gain_sched_parameters, "Correct gain sched model"); - var correct_plantSim = new PlantSimulator(new List { correct_model }); - var inputData = new TimeSeriesDataSet(); - inputData.Add(correct_plantSim.AddExternalSignal(correct_model, SignalType.External_U, (int)INDEX.FIRST), u); - inputData.CreateTimestamps(timeBase_s); - var CorrectisSimulatable = correct_plantSim.Simulate(inputData, out TimeSeriesDataSet CorrectsimData); - SISOTests.CommonAsserts(inputData, CorrectsimData, correct_plantSim); - double[] simY1 = CorrectsimData.GetValues(correct_model.GetID(), SignalType.Output_Y); - unitData.Y_meas = simY1; - - // Act - GainSchedParameters best_params = GainSchedIdentifier.Identify(unitData); - GainSchedModel best_model = new GainSchedModel(best_params, "Best fitting model"); - var best_plantSim = new PlantSimulator(new List { best_model }); - inputData.Add(best_plantSim.AddExternalSignal(best_model, SignalType.External_U, (int)INDEX.FIRST), u); - - var IdentifiedisSimulatable = best_plantSim.Simulate(inputData, out TimeSeriesDataSet IdentifiedsimData); - SISOTests.CommonAsserts(inputData, IdentifiedsimData, best_plantSim); - double[] simY2 = IdentifiedsimData.GetValues(best_model.GetID(), SignalType.Output_Y); - // Number of inputs can be determined from the TimeSeries object, assuming it provides a way to determine this - int numberOfInputs = unitData.U.GetNColumns(); // Example property, replace with actual implementation - - // Assert - int min_number_of_gains = Math.Min(best_params.LinearGainThresholds.Length, correct_gain_sched_parameters.LinearGainThresholds.Length); - for (int k = 0; k < min_number_of_gains; k++) - { - Assert.That(best_params.LinearGainThresholds[k].IsGreaterThanOrEqual(u.First()), - "Linear gain threshold below lower bound (umin) " + ver.ToString()); - Assert.That(best_params.LinearGainThresholds[k].IsLessThanOrEqual(u.Last()), - "Linear gain threshold above upper bound (umax) " + ver.ToString()); - } - - bool doPlot = false; - if (doPlot) - { - Shared.EnablePlots(); - Plot.FromList(new List { - simY1, - simY2, - unitData.U.GetColumn(0) }, - new List { "y1=correct_model", "y1=best_model", "y3=u1" }, - timeBase_s, - "GainSched - Threshold within bounds - " + ver.ToString()); - Shared.DisablePlots(); - } } -/* - [TestCase(1, -0.5)] - [TestCase(2, 2.0)] - [TestCase(3, -1.5)] - [TestCase(4, 3.0)] - [TestCase(5, -3.5)] - [TestCase(6, 4.0)] - [TestCase(7, -4.5)] - public void GainSchedIdentify_AllTimeConstantsArePositive(int ver, double gain_sched_threshold) - { - // 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); - - var gainSchedIdentifier = new GainSchedIdentifier(); - - GainSchedParameters correct_gain_sched_parameters = new GainSchedParameters - { - TimeConstant_s = new double[] { 15, 5 }, - TimeConstantThresholds = new double[] { gain_sched_threshold }, - LinearGains = new List { new double[] { 1 }, new double[] { 3 } }, - LinearGainThresholds = new double[] { gain_sched_threshold }, - TimeDelay_s = 0, - Bias = 0 - }; - GainSchedModel correct_model = new GainSchedModel(correct_gain_sched_parameters, "Correct gain sched model"); - var correct_plantSim = new PlantSimulator(new List { correct_model }); - var inputData = new TimeSeriesDataSet(); - inputData.Add(correct_plantSim.AddExternalSignal(correct_model, SignalType.External_U, (int)INDEX.FIRST), u); - inputData.CreateTimestamps(timeBase_s); - var CorrectisSimulatable = correct_plantSim.Simulate(inputData, out TimeSeriesDataSet CorrectsimData); - SISOTests.CommonAsserts(inputData, CorrectsimData, correct_plantSim); - double[] simY1 = CorrectsimData.GetValues(correct_model.GetID(), SignalType.Output_Y); - unitData.Y_meas = simY1; - - // Act - GainSchedParameters best_params = gainSchedIdentifier.GainSchedIdentify(unitData); - GainSchedModel best_model = new GainSchedModel(best_params, "Best fitting model"); - var best_plantSim = new PlantSimulator(new List { best_model }); - inputData.Add(best_plantSim.AddExternalSignal(best_model, SignalType.External_U, (int)INDEX.FIRST), u); - - var IdentifiedisSimulatable = best_plantSim.Simulate(inputData, out TimeSeriesDataSet IdentifiedsimData); - - SISOTests.CommonAsserts(inputData, IdentifiedsimData, best_plantSim); - - double[] simY2 = IdentifiedsimData.GetValues(best_model.GetID(), SignalType.Output_Y); - - // Number of inputs can be determined from the TimeSeries object, assuming it provides a way to determine this - int numberOfInputs = unitData.U.GetNColumns(); // Example property, replace with actual implementation - - // Assert - int min_number_of_time_constants = best_params.TimeConstant_s.Length; - for (int k = 0; k < min_number_of_time_constants; k++) - { - Assert.That(best_params.TimeConstant_s[k].IsGreaterThanOrEqual((double)0), - "Negative time constant - " + ver.ToString()); - } - - // Shared.EnablePlots(); - Plot.FromList(new List { - simY1, - simY2, - unitData.U.GetColumn(0) }, - new List { "y1=correct_model", "y1=best_model", "y3=u1" }, - timeBase_s, - "GainSched - Positive time constants - " + ver.ToString()); - // Shared.DisablePlots(); - } - */ - - - + } }