Skip to content

Commit

Permalink
- added experimental "dampingratio" to unit model, which when enabled…
Browse files Browse the repository at this point in the history
… gives second order dynamics with overshoot.(not supported by UnitIdentifier yet)
  • Loading branch information
Steinar Elgsæter committed Nov 19, 2024
1 parent 059d265 commit cc5f39d
Show file tree
Hide file tree
Showing 5 changed files with 174 additions and 36 deletions.
35 changes: 12 additions & 23 deletions Dynamic/SimulatableModels/UnitModel.cs
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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)
Expand Down
6 changes: 4 additions & 2 deletions Dynamic/SimulatableModels/UnitParameters.cs
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,11 @@ public class UnitParameters : ModelParametersBaseClass
public double TimeConstant_s { get; set; } = 0;

/// <summary>
/// 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.
/// </summary>
public double TimeConstant2_s { get; set; } = 0;
public double DampingRatio { get; set; } = 0;

/// <summary>
/// The uncertinty of the time constant estimate
Expand Down
34 changes: 24 additions & 10 deletions TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs
Original file line number Diff line number Diff line change
Expand Up @@ -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<ISimulatableModel> { procModel });
new List<ISimulatableModel> { 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<double[]> {
simData.GetValues(procModel.GetID(),SignalType.Output_Y),
simData.GetValues(procModel2.GetID(),SignalType.Output_Y),
inputData.GetValues(procModel.GetID(),SignalType.External_U)},
new List<string> { "y1=y_sim1", "y3=u" },
timeBase_s, TestContext.CurrentContext.Test.Name);
new List<string> { "y1=y_secondorder", "y1=y_firstorder", "y3=u" },
timeBase_s, TestContext.CurrentContext.Test.Name+ v1);
Shared.DisablePlots();
}
CommonAsserts(inputData, simData, plantSim);
Expand Down
2 changes: 1 addition & 1 deletion TimeSeriesAnalysis.csproj
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
<RunAnalyzersDuringLiveAnalysis>False</RunAnalyzersDuringLiveAnalysis>
<RepositoryUrl>https://github.com/equinor/TimeSeriesAnalysis.git</RepositoryUrl>
<PackageReadmeFile>readme.md</PackageReadmeFile>
<Version>1.3.29</Version>
<Version>1.3.30</Version>
<Company>Equinor</Company>
<Authors>Equinor</Authors>
<IncludeSymbols>true</IncludeSymbols>
Expand Down
133 changes: 133 additions & 0 deletions TimeSeriesAnalysis/Filters/SecondOrder.cs
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace TimeSeriesAnalysis
{

///<summary>
/// second-order filter
/// </summary>
public class SecondOrder
{
private double timeBase_s;
private double prevFilteredSignal;
private double prevPrevFilteredSignal;
private int nSignals = 0;
private double nanValue;

/// <summary>
/// Constructor
/// </summary>
/// <param name="TimeBase_s">The time base, the time interval between each time step of the dataset, in seconds</param>
/// <param name="nanValue">value that is to be treated as NaN and ignored</param>
public SecondOrder(double TimeBase_s, double nanValue = -9999)
{
this.timeBase_s = TimeBase_s;
this.nSignals = 0;
this.nanValue = nanValue;
}

/// <summary>
/// Adds a single data point to the filter
/// </summary>
/// <param name="signal">data point</param>
/// <param name="FilterTc_s">filter time constant in seconds</param>
/// <param name="DampingZeta">filter time constant in seconds(also reffered to as zeta )</param>
/// <param name="doReset">usually false, setting to true causes filter to reset to the value of signal</param>
/// <returns></returns>
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;
}

/// <summary>
/// Filter an entire time-series in one command
/// </summary>
/// <param name="signal">the vector of the entire time-series to be filtered</param>
/// <param name="FilterTc_s">filter time constant</param>
/// <param name="order">filter order, either 1 or 2</param>
/// <param name="indicesToIgnore">for these indices the of the signal, the filter should just "freeze" the value(can be null)</param>
/// <returns>a vector of the filtered time-series</returns>
public double[] Filter(double[] signal, double FilterTc_s, int order=1, List<int> 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;
}

}
}

0 comments on commit cc5f39d

Please sign in to comment.