Skip to content

Commit

Permalink
-renaming disturbanceidentifier->disturbancecalculator
Browse files Browse the repository at this point in the history
-fix issues with PlantSimulator.SimulateSingle when using it to back-calculate the disturbance
  • Loading branch information
Steinar Elgsæter committed Dec 17, 2024
1 parent 1e41e50 commit 6516a66
Show file tree
Hide file tree
Showing 5 changed files with 60 additions and 28 deletions.
15 changes: 8 additions & 7 deletions Dynamic/Identification/ClosedLoopUnitIdentifier.cs
Original file line number Diff line number Diff line change
Expand Up @@ -118,12 +118,13 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters

var GSdescription = "";


// ----------------
// step 1: no process model assumed, let disturbance estimator guesstimate a pid-process gain,
// to give a first estimate of the disturbance
var unitModel_step1 = ModelFreeEstimateClosedLoopProcessGain(dataSetRun1, pidParams, pidInputIdx);
var distIdResult1 = DisturbanceIdentifier.CalculateDisturbanceVector(dataSetRun1, unitModel_step1, pidInputIdx, pidParams);
var distIdResult1 = DisturbanceCalculator.CalculateDisturbanceVector(dataSetRun1, unitModel_step1, pidInputIdx, pidParams);
if(distIdResult1.ErrorType != DisturbanceEstimationError.None)
Console.WriteLine("Step 1 Error:" + distIdResult1.ErrorType);
if (doConsoleDebugOut)
ConsoleDebugOut(unitModel_step1, "Step 1");
idUnitModelsList.Add(unitModel_step1);
Expand Down Expand Up @@ -312,7 +313,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters
// approach1:
// "look for the time-constant that gives the disturbance that changes the least"
//
var distIdResult2 = DisturbanceIdentifier.CalculateDisturbanceVector
var distIdResult2 = DisturbanceCalculator.CalculateDisturbanceVector
(dataSetStep3, unitModel, pidInputIdx, pidParams);
var estDisturbances = new List<double[]>();
var distDevs = new List<double>();
Expand All @@ -329,7 +330,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters
var newParams = unitModel.GetModelParameters().CreateCopy();
newParams.TimeConstant_s = candiateTc_s;
var newModel = new UnitModel(newParams);
var distIdResult_Test = DisturbanceIdentifier.CalculateDisturbanceVector
var distIdResult_Test = DisturbanceCalculator.CalculateDisturbanceVector
(dataSetStep3, newModel, pidInputIdx, pidParams);
estDisturbances.Add(distIdResult_Test.d_est);

Expand Down Expand Up @@ -357,7 +358,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters
step3params.TimeConstant_s = candiateTc_s;
var step3Model = new UnitModel(step3params);
idUnitModelsList.Add(step3Model);
var distIdResult_step4 = DisturbanceIdentifier.CalculateDisturbanceVector
var distIdResult_step4 = DisturbanceCalculator.CalculateDisturbanceVector
(dataSetStep3, step3Model, pidInputIdx, pidParams);
idDisturbancesList.Add(distIdResult_step4);
if (doConsoleDebugOut)
Expand Down Expand Up @@ -656,7 +657,7 @@ private static Tuple<UnitModel, double> GlobalSearchLinearPidGain(UnitDataSet da
}
// part 2: actual disturbance modeling
var dataSet_altMISO = new UnitDataSet(dataSet);
var distIdResultAlt_MISO = DisturbanceIdentifier.CalculateDisturbanceVector
var distIdResultAlt_MISO = DisturbanceCalculator.CalculateDisturbanceVector
(dataSet_altMISO, alternativeMISOModel, pidInputIdx, pidParams);
if (false)
{
Expand Down Expand Up @@ -981,7 +982,7 @@ private static Tuple<UnitModel, DisturbanceIdResult> EstimateSISOdisturbanceForP


int pidInputIdx_NewSystem = 0; // pidInputIdx=0 always for the new SISO system:
var candidateModelDisturbance = DisturbanceIdentifier.CalculateDisturbanceVector
var candidateModelDisturbance = DisturbanceCalculator.CalculateDisturbanceVector
(dataSet_SISO, candidateModel, pidInputIdx_NewSystem, pidParams);

var d_est = candidateModelDisturbance.d_est;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,12 @@ namespace TimeSeriesAnalysis.Dynamic
/// <summary>
/// Enum describing reasons that a disturbance may be set by logic to exactly zero
/// </summary>
public enum DisturbanceSetToZeroReason
public enum DisturbanceEstimationError
{
/// <summary>
/// default, when disturbance estimation has not run yet
/// </summary>
NotRunYet=0,
None=0,
/// <summary>
/// Set disturbance to zero beacause simulator did not run
/// </summary>
Expand Down Expand Up @@ -58,7 +58,7 @@ public class DisturbanceIdResult
/// <summary>
/// The reason for setting the disturbance to zero
/// </summary>
public DisturbanceSetToZeroReason zeroReason;
public DisturbanceEstimationError ErrorType;
/// <summary>
///
/// </summary>
Expand Down Expand Up @@ -88,6 +88,7 @@ public class DisturbanceIdResult
public DisturbanceIdResult(UnitDataSet dataSet)
{
N = dataSet.GetNumDataPoints();
ErrorType = DisturbanceEstimationError.None;
SetToZero();
}

Expand All @@ -109,10 +110,9 @@ public void SetToZero()
}

/// <summary>
/// An algorithm that attempts to re-create the additive output disturbance acting on
/// a signal Y while PID-control attempts to counter-act the disturbance by adjusting its manipulated output u.
/// For a given process model and dataset, calculate the disturbance vector d by subtracting y_proc from y_meas
/// </summary>
public class DisturbanceIdentifier
public class DisturbanceCalculator
{

/// <summary>
Expand Down Expand Up @@ -293,7 +293,7 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat

if (y_proc == null)
{
result.zeroReason = DisturbanceSetToZeroReason.UnitSimulatorUnableToRun;
result.ErrorType = DisturbanceEstimationError.UnitSimulatorUnableToRun;
return result;
}

Expand Down
52 changes: 41 additions & 11 deletions Dynamic/PlantSimulator/PlantSimulator.cs
Original file line number Diff line number Diff line change
Expand Up @@ -543,18 +543,24 @@ public static (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl

/// <summary>
/// Simulate single model based on a unit data set
///
/// This is a convenience function that creates a TimeSeriesDataSet, sets default names in the model and dataset that match based on unitDataset
/// The output is returned directly.
///
/// Optionally, the result can be written to y_meas or y_sim in unitdata.
/// </summary>
/// <param name="unitData">contains a unit data set that must have U filled, Y_sim will be written here</param>
/// <param name="model">model to simulate</param>
/// <param name="writeToYmeas">if set to true, the simulated result is written to unitData.Y_meas instead of Y_sim</param>
/// <param name="noiseAmplitude">if writing to Ymeas, it is possible to add noise of the given amplitude to signal</param>
/// <param name="addSimToUnitData">if true, the Y_sim of unitData has the simulation result written two i</param>
/// <param name="seedNr">the seed value of the noise to be added</param>
/// <returns>a tuple, first aa true if able to simulate, otherwise false, second is the simulated time-series</returns>
static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatableModel model,bool writeToYmeas= false,
/// <returns>a tuple, first aa true if able to simulate, otherwise false, second is the simulated time-series "y_proc" without any additive </returns>
static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatableModel model, bool writeToYmeas= false,
double noiseAmplitude=0,
bool addSimToUnitData=false, int seedNr=123)
{
const string defaultOutputName = "output";
var inputData = new TimeSeriesDataSet();
var singleModelName = "SimulateSingle";
var modelCopy = model.Clone(singleModelName);
Expand All @@ -573,18 +579,42 @@ static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl
inputData.Add(uName, unitData.U.GetColumn(colIdx));
uNames.Add(uName);
}
modelCopy.SetInputIDs(uNames.ToArray());
if (modelCopy.GetOutputID() == null)
modelCopy.SetOutputID("output");

var sim = new PlantSimulator(new List<ISimulatableModel> { modelCopy });
var isOk = sim.SimulateSingle(inputData, singleModelName, false, out var simData);
modelCopy.SetInputIDs(uNames.ToArray());
{
inputData.Add(defaultOutputName, unitData.Y_meas);
modelCopy.SetOutputID(defaultOutputName);
}

// var isOk = PlantSimulator.SimulateSingle(inputData, modelCopy, out var simData);
var sim = new PlantSimulator(new List<ISimulatableModel> { modelCopy });
var isOk = sim.SimulateSingle(inputData, singleModelName, false, out var simData);

if(!isOk)
return (false, null);
double[] y_sim = simData.GetValues(singleModelName, SignalType.Output_Y);

double[] y_proc = null;
double[] y_sim = null;

y_sim = simData.GetValues(defaultOutputName);

if (simData.ContainsSignal(singleModelName))
{
y_proc = simData.GetValues(singleModelName);
}
else
{
y_proc = y_sim;
}


// if the input included a "additive distubance" signal,there will be a internal process output
/* double[] y_proc = simData.GetValues(singleModelName);
double[] y_sim = y_proc;
if (y_proc == null)
{
y_proc = simData.GetValues(defaultOutputName);
y_sim = simData.GetValues(singleModelName, SignalType.Output_Y);
}*/
if (noiseAmplitude > 0)
{
// use a specific seed here, to avoid potential issues with "random unit tests" and not-repeatable
Expand All @@ -607,7 +637,7 @@ static private (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatabl
unitData.Y_sim = y_sim;
}
}
return (isOk, y_sim);
return (isOk, y_proc);
}


Expand Down
5 changes: 3 additions & 2 deletions TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,8 @@ public void Static_StepDisturbance_EstimatesOk(double stepAmplitude)
[TestCase(5)]
public void Static_SinusDisturbance_EstimatesOk(double stepAmplitude)
{
var trueDisturbance = TimeSeriesCreator.Sinus(stepAmplitude, timeBase_s*15, timeBase_s,N );
var sinusPeriod = timeBase_s * 15;
var trueDisturbance = TimeSeriesCreator.Sinus(stepAmplitude, sinusPeriod, timeBase_s,N );
GenericDisturbanceTest(new UnitModel(staticModelParameters, "StaticProcess"), trueDisturbance);
}

Expand Down Expand Up @@ -149,7 +150,7 @@ public void GenericDisturbanceTest (UnitModel processModel, double[] trueDistur
Assert.IsTrue(isOk);
Assert.IsTrue(simData.ContainsSignal(processModel.GetID()),"simulated dataset should include internal process model output (pre-disturbance)");
var pidDataSet = plantSim.GetUnitDataSetForPID(inputData.Combine(simData), pidModel1);
var result = DisturbanceIdentifier.CalculateDisturbanceVector(pidDataSet, processModel);
var result = DisturbanceCalculator.CalculateDisturbanceVector(pidDataSet, processModel);


if (doAssertResult)
Expand Down
2 changes: 1 addition & 1 deletion docs/articles/processsimulator.md
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ Each model has a number of inputs that travel through the model, each with an id
In addition, models support adding signals directly to the output. This is a feature intended for modeling disturbances.
The IDs of such signals are referred to as *additive input IDs*.

### Estimating disturbances and simulating their effects
### Estimating disturbances and simulating their effects using PlantSimulator.SimulateSingle()

If the ``inputData`` given to the PlantSimulator includes measured pid-outputs ``u`` and process outputs ``y`` of feedback loops,
then PlantSimulator uses that
Expand Down

0 comments on commit 6516a66

Please sign in to comment.