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;
+ }
+
+ }
+}