From cc5f39d9466c9ece6b5facc635d9621c071a3176 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steinar=20Elgs=C3=A6ter?= Date: Tue, 19 Nov 2024 15:00:43 +0100 Subject: [PATCH] - added experimental "dampingratio" to unit model, which when enabled gives second order dynamics with overshoot.(not supported by UnitIdentifier yet) --- Dynamic/SimulatableModels/UnitModel.cs | 35 ++--- Dynamic/SimulatableModels/UnitParameters.cs | 6 +- .../Tests/PlantSimulatorSISOTests.cs | 34 +++-- TimeSeriesAnalysis.csproj | 2 +- TimeSeriesAnalysis/Filters/SecondOrder.cs | 133 ++++++++++++++++++ 5 files changed, 174 insertions(+), 36 deletions(-) create mode 100644 TimeSeriesAnalysis/Filters/SecondOrder.cs diff --git a/Dynamic/SimulatableModels/UnitModel.cs b/Dynamic/SimulatableModels/UnitModel.cs index 830602bf..c7c01cdd 100644 --- a/Dynamic/SimulatableModels/UnitModel.cs +++ b/Dynamic/SimulatableModels/UnitModel.cs @@ -53,10 +53,10 @@ public class UnitModel : ModelBaseClass, ISimulatableModel private LowPass lowPass1order; // second order dynamics : - private LowPass lowPass2order1; - private LowPass lowPass2order2; - private TimeDelay delayObj2order; + private SecondOrder secondOrder; + + // time-delay private TimeDelay delayObj; @@ -505,30 +505,19 @@ public double[] Iterate(double[] inputs, double timeBase_s,double badValueIndica // nb! if first iteration, start model at steady-state double x_dynamic = x_ss; - if (modelParameters.TimeConstant_s >= 0) + if (modelParameters.TimeConstant_s >= 0 && modelParameters.DampingRatio == 0) { x_dynamic = lowPass1order.Filter(x_ss, modelParameters.TimeConstant_s, 1, isFirstIteration); - - // second order time constant(work in progress) - /* if (modelParameters.TimeConstant2_s > 0) + } + else if (modelParameters.TimeConstant_s >= 0 && modelParameters.DampingRatio > 0) + { + if (this.secondOrder == null) { - if (this.lowPass2order1 == null) - { - this.lowPass2order1 = new LowPass(timeBase_s); - } - if (this.delayObj2order == null) - { - this.delayObj2order = new TimeDelay(timeBase_s, timeBase_s); - } - - double secondOrderFilt = lowPass2order1.Filter(x_ss, modelParameters.TimeConstant2_s, 1, isFirstIteration); - double delayed_secondOrderFilt = delayObj2order.Delay(secondOrderFilt); - - double x_secondOrderDyn = (secondOrderFilt- delayed_secondOrderFilt)/timeBase_s; - x_dynamic -= x_secondOrderDyn; - }*/ - + this.secondOrder = new SecondOrder(timeBase_s); + } + x_dynamic = secondOrder.Filter(x_ss, modelParameters.TimeConstant_s, modelParameters.DampingRatio, isFirstIteration); } + isFirstIteration = false; double y = 0; if (modelParameters.TimeDelay_s <= 0) diff --git a/Dynamic/SimulatableModels/UnitParameters.cs b/Dynamic/SimulatableModels/UnitParameters.cs index 647bbf55..bb7a3680 100644 --- a/Dynamic/SimulatableModels/UnitParameters.cs +++ b/Dynamic/SimulatableModels/UnitParameters.cs @@ -33,9 +33,11 @@ public class UnitParameters : ModelParametersBaseClass public double TimeConstant_s { get; set; } = 0; /// - /// A time constant in seconds, the second order time constant. Leave to zero to turn off second order dynamics. + /// Damping (second-order) values between ~0.3-0.99 will cause step response with a single visibl overshoot. ) + /// Set to zero to disable damping. + /// As values less than 0.3 approach zero, the step response will become more and more oscillatory. /// - public double TimeConstant2_s { get; set; } = 0; + public double DampingRatio { get; set; } = 0; /// /// The uncertinty of the time constant estimate diff --git a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs b/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs index 03118ee0..96833590 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs @@ -115,35 +115,49 @@ public void SimulateSingle_InitsRunsAndConverges() [TestCase,Explicit] public void SimulateSingle_SecondOrderSystem() { - var mp = new UnitParameters + string v1 = "6"; + int N = 5000; + int Nstep = 50; + + var modelParams2ndOrder = new UnitParameters { - TimeConstant_s = 10, - TimeConstant2_s = 2, + TimeConstant_s = 50, + DampingRatio = 0.1,// seems that if dampingRatio is higher than 1, it is redundant!(coudl be described ) + LinearGains = new double[] { 1 }, + TimeDelay_s = 0, + Bias = 5 + }; + var modelParams1stOrder = new UnitParameters + { + TimeConstant_s = 50, + DampingRatio = 0, LinearGains = new double[] { 1 }, TimeDelay_s = 0, Bias = 5 }; - var procModel = new UnitModel(mp,"second order system"); + var procModel = new UnitModel(modelParams2ndOrder,"second order system"); + var procModel2 = new UnitModel(modelParams1stOrder, "first order system"); var plantSim = new PlantSimulator( - new List { procModel }); + new List { procModel,procModel2 }); var inputData = new TimeSeriesDataSet(); - inputData.Add(plantSim.AddExternalSignal(procModel, SignalType.External_U), TimeSeriesCreator.Step(N / 4, N, 50, 55)); + inputData.Add(plantSim.AddExternalSignal(procModel, SignalType.External_U), TimeSeriesCreator.Step(Nstep, N, 50, 55)); + inputData.Add(plantSim.AddExternalSignal(procModel2, SignalType.External_U), TimeSeriesCreator.Step(Nstep, N, 50, 55)); inputData.CreateTimestamps(timeBase_s); var isOk = plantSim.Simulate(inputData,out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - - if (true) + if (false) { Shared.EnablePlots(); Plot.FromList(new List { simData.GetValues(procModel.GetID(),SignalType.Output_Y), + simData.GetValues(procModel2.GetID(),SignalType.Output_Y), inputData.GetValues(procModel.GetID(),SignalType.External_U)}, - new List { "y1=y_sim1", "y3=u" }, - timeBase_s, TestContext.CurrentContext.Test.Name); + new List { "y1=y_secondorder", "y1=y_firstorder", "y3=u" }, + timeBase_s, TestContext.CurrentContext.Test.Name+ v1); Shared.DisablePlots(); } CommonAsserts(inputData, simData, plantSim); diff --git a/TimeSeriesAnalysis.csproj b/TimeSeriesAnalysis.csproj index ab59e07c..eb2be044 100644 --- a/TimeSeriesAnalysis.csproj +++ b/TimeSeriesAnalysis.csproj @@ -14,7 +14,7 @@ False https://github.com/equinor/TimeSeriesAnalysis.git readme.md - 1.3.29 + 1.3.30 Equinor Equinor true diff --git a/TimeSeriesAnalysis/Filters/SecondOrder.cs b/TimeSeriesAnalysis/Filters/SecondOrder.cs new file mode 100644 index 00000000..c00397d3 --- /dev/null +++ b/TimeSeriesAnalysis/Filters/SecondOrder.cs @@ -0,0 +1,133 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; + +namespace TimeSeriesAnalysis +{ + + /// + /// second-order filter + /// + public class SecondOrder + { + private double timeBase_s; + private double prevFilteredSignal; + private double prevPrevFilteredSignal; + private int nSignals = 0; + private double nanValue; + + /// + /// Constructor + /// + /// The time base, the time interval between each time step of the dataset, in seconds + /// value that is to be treated as NaN and ignored + public SecondOrder(double TimeBase_s, double nanValue = -9999) + { + this.timeBase_s = TimeBase_s; + this.nSignals = 0; + this.nanValue = nanValue; + } + + /// + /// Adds a single data point to the filter + /// + /// data point + /// filter time constant in seconds + /// filter time constant in seconds(also reffered to as zeta ) + /// usually false, setting to true causes filter to reset to the value of signal + /// + public double Filter(double signal, double FilterTc_s, double DampingZeta=0, bool doReset = false) + { + // DampingZeta = 0; + + double startUpSignals = 2; + if (nSignals < startUpSignals) + { + nSignals++; + this.prevPrevFilteredSignal = signal;// (*Force filter to steady - state *) + this.prevFilteredSignal = signal; // (*Force filter to steady - state *) + } + //double a; + double filteredSignal= signal; + double omega_n = 1 / FilterTc_s;//todo: is this correct? + + double divisor = (1 / timeBase_s + 2 * DampingZeta * omega_n / timeBase_s + Math.Pow(omega_n,2)); + if (divisor <= 0) + return signal; + + double a_1 = (2 + 2 * DampingZeta * omega_n) / timeBase_s / divisor; + double a_2 = -1 / timeBase_s / divisor; + double b = 1-a_1-a_2 ; + + if (Double.IsNaN(this.prevFilteredSignal) || Double.IsInfinity(this.prevFilteredSignal) + || prevFilteredSignal == nanValue || Double.IsNaN(this.prevPrevFilteredSignal) + || Double.IsInfinity(this.prevPrevFilteredSignal) + || prevPrevFilteredSignal == nanValue) + filteredSignal = signal; + else + { + filteredSignal = a_1 * this.prevFilteredSignal + a_2 * this.prevPrevFilteredSignal + b* signal; + } + this.prevPrevFilteredSignal = this.prevFilteredSignal; + this.prevFilteredSignal = filteredSignal; + // if (nSignals <= startUpSignals-1) + // doReset = true; + if (doReset) + { + this.prevPrevFilteredSignal = signal;// (*Force filter to steady - state *) + this.prevFilteredSignal = signal; // (*Force filter to steady - state *) + filteredSignal = signal; + } + return filteredSignal; + } + + /// + /// Filter an entire time-series in one command + /// + /// the vector of the entire time-series to be filtered + /// filter time constant + /// filter order, either 1 or 2 + /// for these indices the of the signal, the filter should just "freeze" the value(can be null) + /// a vector of the filtered time-series + public double[] Filter(double[] signal, double FilterTc_s, int order=1, List indicesToIgnore= null) + { + double[] outSig = new double[signal.Length]; + + if (indicesToIgnore == null) + { + for (int i = 0; i < signal.Count(); i++) + { + outSig[i] = this.Filter(signal[i], FilterTc_s, order, false); + } + } + else + { + int lastGoodIndex = 0; + for (int i = 0; i < signal.Count(); i++) + { + if (indicesToIgnore.Contains(i)) + { + if (i == 0) + { + int i_forward = 1; + while (i_forward < signal.Count()-1 && indicesToIgnore.Contains(i_forward)) + { + i_forward++; + } + lastGoodIndex = i_forward; + } + } + else + { + lastGoodIndex = i; + } + outSig[i] = this.Filter(signal[lastGoodIndex], FilterTc_s, order, false); + } + } + return outSig; + } + + } +}