diff --git a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs index 0fd83ea7..61ec9d9a 100644 --- a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs +++ b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs @@ -138,7 +138,15 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters idDisturbancesList.Add(distIdResult1); idUnitModelsList.Add(unitModel_step1); if (doConsoleDebugOut) - ConsoleDebugOut(unitModel_step1, "Step 1,MISO"); + { + string otherGains = ""; + for(int i=0;i < unitModel_step1.GetModelParameters().LinearGains.Count();i++) + { + if (i!= pidInputIdx) + otherGains+= " other Gains: "+ unitModel_step1.GetModelParameters().LinearGains.ElementAt(i).ToString("F2"); + } + ConsoleDebugOut(unitModel_step1, "Step 1,MISO", otherGains); + } } // EstimateDisturbanceLF(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); @@ -213,8 +221,6 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters } } - - // ---------------- // - issue is that after run 1 modelled output does not match measurement. // - the reason that we want this run is that after run1 the time constant and @@ -260,7 +266,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters DateTime[] dateTimes = null; bool doDebugPlot = true; - (var plantSim, var inputData) = PlantSimulator.CreateFeedbackLoopWithEstimatedDisturbance(dataSetStep3, pidModel, unitModel, pidInputIdx); + (var plantSim, var inputData) = PlantSimulatorHelper.CreateFeedbackLoopWithEstimatedDisturbance(dataSetStep3, pidModel, unitModel, pidInputIdx); if (doDebugPlot) { y_simList.Add(inputData.GetValues(unitModel.ID, SignalType.Output_Y)); @@ -459,7 +465,7 @@ private static double EstimateDisturbanceLF(UnitDataSet dataSet, UnitModel unitM { unitParams.LinearGains = new double[] { Kp }; umInternal.SetModelParameters(unitParams); - (var isOk, var y_proc) = PlantSimulator.SimulateSingle(dataSet, umInternal); + (var isOk, var y_proc) = PlantSimulatorHelper.SimulateSingle(dataSet, umInternal); var d_LF = vec.Multiply(vec.Subtract(y_proc, y_proc[0]), -1); var d_est1 = vec.Add(d_HF, d_LF); var d_est2 = vec.Subtract(dataSet.Y_meas, y_proc); @@ -1018,7 +1024,7 @@ public static bool ClosedLoopSim(UnitDataSet unitData, UnitParameters modelParam //TODO: this does not ignore bad datapoints? var pidModel = new PidModel(pidParams,"PidModel"); var unitModel = new UnitModel(modelParams, "ProcModel"); - (var plantSim, var inputDataSet) = PlantSimulator.CreateFeedbackLoopWithEstimatedDisturbance(unitData, pidModel, + (var plantSim, var inputDataSet) = PlantSimulatorHelper.CreateFeedbackLoopWithEstimatedDisturbance(unitData, pidModel, unitModel, pidInputIdx); var isOk = plantSim.Simulate(inputDataSet, out var simData); diff --git a/Dynamic/Identification/DisturbanceCalculator.cs b/Dynamic/Identification/DisturbanceCalculator.cs index 986f36e7..865f20ba 100644 --- a/Dynamic/Identification/DisturbanceCalculator.cs +++ b/Dynamic/Identification/DisturbanceCalculator.cs @@ -17,6 +17,7 @@ using TimeSeriesAnalysis; using TimeSeriesAnalysis.Utility; +using Newtonsoft.Json.Linq; namespace TimeSeriesAnalysis.Dynamic @@ -72,14 +73,6 @@ public class DisturbanceIdResult /// (For debugging:) an estimated process gain /// public double estPidProcessGain; - /// - /// (for debugging) the high-frequency component of the disturbance that is determined by observing changes in process variable y - /// - public double[] d_HF; - /// - /// (for debugging) the low-frequency component of the disturbance that is determined by observing changes in manipulated variable u - /// - public double[] d_LF; /// /// Constuctor @@ -100,13 +93,9 @@ public void SetToZero() d_est = Vec.Fill(0, N); estPidProcessGain = 0; isAllZero = true; - d_HF = Vec.Fill(0, N); - d_LF = Vec.Fill(0, N); adjustedUnitDataSet = null; - } - } /// @@ -158,12 +147,14 @@ private static UnitDataSet RemoveSetpointAndOtherInputChangeEffectsFromDataSet(U // BEGIN "no_dist" process simulation = // a simulation of the process that does not include any real Y_meas or u_pid, thus no effects of // disturbances are visible in this simulation - - + + var unitModelCopy = (UnitModel)unitModel.Clone("RemoveSetpointAndOtherInputChangeEffectsFromDataSet");// make copy that has no additive output signals(disturbance) + unitModelCopy.additiveInputIDs = null; + var processSim_noDist = new PlantSimulator( - new List { pidModel1, unitModel }); - processSim_noDist.ConnectModels(unitModel, pidModel1); - processSim_noDist.ConnectModels(pidModel1, unitModel,pidInputIdx); + new List { pidModel1, unitModelCopy }); + processSim_noDist.ConnectModels(unitModelCopy, pidModel1); + processSim_noDist.ConnectModels(pidModel1, unitModelCopy, pidInputIdx); var inputData_noDist = new TimeSeriesDataSet(); if (unitDataSet.U.GetNColumns()>1) @@ -172,7 +163,7 @@ private static UnitDataSet RemoveSetpointAndOtherInputChangeEffectsFromDataSet(U { if (curColIdx == pidInputIdx) continue; - inputData_noDist.Add(processSim_noDist.AddExternalSignal(unitModel, SignalType.External_U, curColIdx), + inputData_noDist.Add(processSim_noDist.AddExternalSignal(unitModelCopy, SignalType.External_U, curColIdx), unitDataSet.U.GetColumn(curColIdx)); } } @@ -203,7 +194,8 @@ private static UnitDataSet RemoveSetpointAndOtherInputChangeEffectsFromDataSet(U // create a new Y_meas that excludes the influence of any disturbance using "no_Dist" simulation // this is used to find d_HF - var procOutputY = simData_noDist.GetValues(unitModel.GetID(), SignalType.Output_Y); + //var procOutputY = simData_noDist.GetValues(unitModel.GetID(), SignalType.Output_Y); + var procOutputY = simData_noDist.GetValues(unitModelCopy.GetID()); var deltaProcOutputY = vec.Subtract(procOutputY, procOutputY[idxFirstGoodValue]); unitDataSet_adjusted.Y_meas = vec.Subtract(unitDataSet.Y_meas, deltaProcOutputY); @@ -289,7 +281,7 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat // non-disturbance related changes in the dataset producing "unitDataSet_adjusted" var unitDataSet_adjusted = RemoveSetpointAndOtherInputChangeEffectsFromDataSet(unitDataSet, unitModel, pidInputIdx, pidParams); unitDataSet_adjusted.D = null; - (bool isOk, double[] y_proc) = PlantSimulator.SimulateSingle(unitDataSet_adjusted, unitModel); + (bool isOk, double[] y_proc) = PlantSimulatorHelper.SimulateSingle(unitDataSet_adjusted, unitModel); if (y_proc == null) { @@ -311,13 +303,35 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat } //TODO: can these be removed? - double[] d_LF = vec.Multiply(vec.Subtract(y_proc, y_proc[indexOfFirstGoodValue]), -1); - double[] d_HF = vec.Subtract(unitDataSet_adjusted.Y_meas, unitDataSet_adjusted.Y_setpoint); + //double[] d_LF = vec.Multiply(vec.Subtract(y_proc, y_proc[indexOfFirstGoodValue]), -1); + //double[] d_HF = vec.Subtract(unitDataSet_adjusted.Y_meas, unitDataSet_adjusted.Y_setpoint); + + // old: d[k] = y_meas[k] -y_proc[k] + //double[] d_est = vec.Subtract(unitDataSet_adjusted.Y_meas, y_proc); + + bool IsNaN(double value) + { + if (double.IsNaN(value) || value == unitDataSet_adjusted.BadDataID) + return true; + else + return false; + } + + double[] d_est = new double[unitDataSet_adjusted.Y_meas.Length]; + + for (int i = 1; i < d_est.Length; i++) + { + if (IsNaN(unitDataSet_adjusted.Y_meas[i]) || IsNaN(y_proc[i])) + d_est[i] = double.NaN; + else + d_est[i] = unitDataSet_adjusted.Y_meas[i] - y_proc[i-1]; + } + d_est[0] = unitDataSet_adjusted.Y_meas[0] - y_proc[0]; + - double[] d_est = vec.Subtract(unitDataSet_adjusted.Y_meas, y_proc); result.d_est = d_est; - result.d_LF = d_LF; - result.d_HF = d_HF; + //result.d_LF = d_LF; + // result.d_HF = d_HF; result.estPidProcessGain = unitModel.GetModelParameters().LinearGains.ElementAt(pidInputIdx); result.adjustedUnitDataSet = unitDataSet_adjusted; diff --git a/Dynamic/Identification/GainSchedIdentifier.cs b/Dynamic/Identification/GainSchedIdentifier.cs index e19eebc9..15d91e9d 100644 --- a/Dynamic/Identification/GainSchedIdentifier.cs +++ b/Dynamic/Identification/GainSchedIdentifier.cs @@ -388,7 +388,7 @@ static private (GainSchedParameters, int) ChooseBestModelFromFittingInfo( if (simulateAndAddToYsimInDataSet) { var bestModel = new GainSchedModel(BestGainSchedParams); - PlantSimulator.SimulateSingleToYsim(dataSet, bestModel); + PlantSimulatorHelper.SimulateSingleToYsim(dataSet, bestModel); } return (BestGainSchedParams, bestModelIdx); } @@ -512,7 +512,7 @@ private static bool DetermineOperatingPointAndSimulate(ref GainSchedParameters g } gsParams.MoveOperatingPointUWithoutChangingModel(desiredOpU); - (var isOk, var y_sim) = PlantSimulator.SimulateSingle(dataSet, gsIdentModel); + (var isOk, var y_sim) = PlantSimulatorHelper.SimulateSingle(dataSet, gsIdentModel); if (isOk) { @@ -523,7 +523,7 @@ private static bool DetermineOperatingPointAndSimulate(ref GainSchedParameters g if (estBias.HasValue) { gsParams.IncreaseOperatingPointY(estBias.Value); - (var isOk2, var y_sim2) = PlantSimulator.SimulateSingle(dataSet, gsIdentModel); + (var isOk2, var y_sim2) = PlantSimulatorHelper.SimulateSingle(dataSet, gsIdentModel); dataSet.Y_sim = y_sim2; if (gsParams.Fitting == null) gsParams.Fitting = new FittingInfo(); @@ -587,7 +587,7 @@ private static void EstimateTimeDelay(ref GainSchedParameters gsParams, ref Unit copiedGsParams.TimeConstant_s = vec.Subtract(gsParams.TimeConstant_s, timedelay_s); var gsIdentModel = new GainSchedModel(copiedGsParams, "ident_model"); - (var isOk, var y_sim) = PlantSimulator.SimulateSingle(dataSet, gsIdentModel); + (var isOk, var y_sim) = PlantSimulatorHelper.SimulateSingle(dataSet, gsIdentModel); if (isOk) { diff --git a/Dynamic/Identification/UnitIdentifier.cs b/Dynamic/Identification/UnitIdentifier.cs index c78d353f..73c26803 100644 --- a/Dynamic/Identification/UnitIdentifier.cs +++ b/Dynamic/Identification/UnitIdentifier.cs @@ -118,7 +118,7 @@ public static UnitModel IdentifyLinearDiff(ref UnitDataSet dataSet, FittingSpecs if (model.modelParameters.Fitting.WasAbleToIdentify) { - PlantSimulator.SimulateSingleToYsim(dataSet, model); + PlantSimulatorHelper.SimulateSingleToYsim(dataSet, model); model.SetFittedDataSet(dataSet); } return model; @@ -449,7 +449,7 @@ private static UnitModel Identify_Internal(ref UnitDataSet dataSet, FittingSpecs // simulate if (modelParameters.Fitting.WasAbleToIdentify) { - PlantSimulator.SimulateSingleToYsim(dataSet, model); + PlantSimulatorHelper.SimulateSingleToYsim(dataSet, model); model.SetFittedDataSet(dataSet); } return model; diff --git a/Dynamic/PlantSimulator/Comment.cs b/Dynamic/PlantSimulator/Comment.cs new file mode 100644 index 00000000..b6e67b3c --- /dev/null +++ b/Dynamic/PlantSimulator/Comment.cs @@ -0,0 +1,44 @@ +using System; +using System.Collections.Generic; +using System.Text; + +namespace TimeSeriesAnalysis.Dynamic +{ + /// + /// Class that holds comments added to models. + /// + public class Comment + { + /// + /// Author of comment. + /// + public string author; + /// + /// Date of comment. + /// + public DateTime date; + /// + /// Comment string + /// + public string comment; + /// + /// Plant score, intended to hold manully set values indicating specific statuses of the model. + /// + public double plantScore; + + /// + /// Comment constructor. + /// + /// + /// + /// + /// + public Comment(string author, DateTime date, string comment, double plantScore = 0) + { + this.author = author; + this.date = date; + this.comment = comment; + this.plantScore = plantScore; + } + } +} diff --git a/Dynamic/PlantSimulator/PlantSimulator.cs b/Dynamic/PlantSimulator/PlantSimulator.cs index c62eebfe..5aafc58b 100644 --- a/Dynamic/PlantSimulator/PlantSimulator.cs +++ b/Dynamic/PlantSimulator/PlantSimulator.cs @@ -22,43 +22,7 @@ namespace TimeSeriesAnalysis.Dynamic { - /// - /// Class that holds comments added to models. - /// - public class Comment - { - /// - /// Author of comment. - /// - public string author; - /// - /// Date of comment. - /// - public DateTime date; - /// - /// Comment string - /// - public string comment; - /// - /// Plant score, intended to hold manully set values indicating specific statuses of the model. - /// - public double plantScore; - /// - /// Comment constructor. - /// - /// - /// - /// - /// - public Comment(string author, DateTime date, string comment, double plantScore =0) - { - this.author = author; - this.date = date; - this.comment = comment; - this.plantScore= plantScore; - } - } /// /// Simulates larger "plant-models" that is built up of connected sub-models @@ -318,88 +282,35 @@ public string ConnectModels(ISimulatableModel upstreamModel, ISimulatableModel d } /// - /// Create a PlantSimulator and TimeSeriesDataSet from a UnitDataSet, PidModel and UnitModel to do closed-loop simulations - /// - /// The feedback loop has NO disturbance signal added, but this can be added to the returned PlantSimulator as needed. - /// + /// Get a list of all the outputs that are pid-controlled and thus need different treatment in the simualor. + /// + /// NOte that in the case that a PIDmodel is simulated in isolation using for example SimulateSingle, the /// - /// - /// - /// - /// - /// a simulator object and a dataset object that is ready to be simulated with Simulate() - public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoopNoDisturbance(UnitDataSet unitDataSet, PidModel pidModel, - UnitModel unitModel, int pidInputIdx=0) + /// a dictionary of the output signal ids and the modelIDs used to generate them. ModelIDs can be null + private Dictionary DeterminePidControlledOutputs() { - var plantSim = new PlantSimulator( - new List { pidModel, unitModel }); - var signalId1 = plantSim.ConnectModels(unitModel, pidModel); - var signalId2 = plantSim.ConnectModels(pidModel, unitModel, pidInputIdx); - - var inputData = new TimeSeriesDataSet(); - inputData.Add(signalId1, (double[])unitDataSet.Y_meas.Clone()); + var ret = new Dictionary(); - for (int curColIdx = 0; curColIdx < unitDataSet.U.GetNColumns(); curColIdx++) + foreach (var model in modelDict) { - if (curColIdx == pidInputIdx) + if (model.Value.GetProcessModelType() == ModelType.PID) { - inputData.Add(signalId2, (double[])unitDataSet.U.GetColumn(pidInputIdx).Clone()); - } - else - { - inputData.Add(plantSim.AddExternalSignal(unitModel, SignalType.External_U, curColIdx), - (double[])unitDataSet.U.GetColumn(curColIdx).Clone()); + var pidControlledOutputSignalID = model.Value.GetModelInputIDs().ElementAt((int)PidModelInputsIdx.Y_meas); + string modelThatCreatesOutputID = null; + foreach (var otherModel in modelDict) + { + if (otherModel.Value.GetOutputID() == pidControlledOutputSignalID) + { + modelThatCreatesOutputID = otherModel.Value.GetID(); + } + } + // if(modelThatCreatesOutputID != null) + ret.Add(model.Value.GetModelInputIDs().ElementAt((int)PidModelInputsIdx.Y_meas), modelThatCreatesOutputID); } } - inputData.Add(plantSim.AddExternalSignal(pidModel, SignalType.Setpoint_Yset), (double[])unitDataSet.Y_setpoint.Clone()); - inputData.CreateTimestamps(unitDataSet.GetTimeBase()); - inputData.SetIndicesToIgnore(unitDataSet.IndicesToIgnore); - return (plantSim, inputData); - } - - /// - /// Create a feedback loop, where the process model has an additive disturbance that is to be estimated. - /// - /// - /// - /// - /// - /// a simulator object and a dataset object that is ready to be simulated with Simulate() - public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoopWithEstimatedDisturbance(UnitDataSet unitDataSet, PidModel pidModel, - UnitModel unitModel, int pidInputIdx = 0) - { - // vital that signal follows naming convention, otherwise it will not be estimated, but should be provided. - unitModel.AddSignalToOutput(SignalNamer.EstDisturbance(unitModel)); - (var sim, var data) = CreateFeedbackLoopNoDisturbance(unitDataSet,pidModel, unitModel, pidInputIdx); - return (sim,data); - } - - - /// - /// Returns a unit data set for a given UnitModel. - /// - /// - /// - /// - public UnitDataSet GetUnitDataSetForProcess(TimeSeriesDataSet inputData, UnitModel unitModel) - { - UnitDataSet dataset = new UnitDataSet(); - dataset.U = new double[inputData.GetLength().Value, 1]; - - dataset.Times = inputData.GetTimeStamps(); - var inputIDs = unitModel.GetModelInputIDs(); - var outputID = unitModel.GetOutputID(); - dataset.Y_meas = inputData.GetValues(outputID); - for (int inputIDidx = 0; inputIDidx < inputIDs.Length; inputIDidx++) - { - var inputID = inputIDs[inputIDidx]; - var curCol = inputData.GetValues(inputID); - dataset.U.WriteColumn(inputIDidx, curCol); - } - return dataset; + return ret; } - /// /// Returns a "unitDataSet" for the given pidModel in the plant. /// This function only works when the unit model connected to the pidModel only has a single input. @@ -539,18 +450,18 @@ double[] GetValuesFromEitherDatasetInternal( int timeIndexInternal) if (model.GetProcessModelType() == ModelType.PID && timeIndex > 0) { int lookBackIndex = 1; - double[] lookBackValues = GetValuesFromEitherDatasetInternal( timeIndex - lookBackIndex); ; + // double[] lookBackValues = GetValuesFromEitherDatasetInternal( timeIndex - lookBackIndex); ; double[] currentValues = GetValuesFromEitherDatasetInternal( timeIndex); // "use values from current data point when available, but fall back on using values from the previous sample if need be" // for instance, always use the most current setpoint value, but if no disturbance vector is given, then use the y_proc simulated from the last iteration. double[] retValues = new double[currentValues.Length]; - retValues = lookBackValues; + retValues = currentValues; // adding in the below code seems to remove the issue with there being a one sample wait time before the effect of a setpoint // is seen on the output, but causes there to be small deviation between what the PlantSimulator.SimulateSingle and PlantSimulator.Simulate // seem to return for a PID-loop in the test BasicPID_CompareSimulateAndSimulateSingle_MustGiveSameResultForDisturbanceEstToWork - for (int i = 0; i < currentValues.Length; i++) + /*for (int i = 0; i < currentValues.Length; i++) { if (Double.IsNaN(currentValues[i])) { @@ -560,7 +471,7 @@ double[] GetValuesFromEitherDatasetInternal( int timeIndexInternal) { retValues[i] = currentValues[i]; } - } + }*/ return retValues; } else @@ -590,138 +501,12 @@ public bool SimulateSingleWithoutAdditive(TimeSeriesDataSet inputData, string si /// /// /// - public bool SimulateSingle(TimeSeriesDataSet inputData, string singleModelName, out TimeSeriesDataSet simData) + /* public bool SimulateSingle(TimeSeriesDataSet inputData, string singleModelName, out TimeSeriesDataSet simData) { return SimulateSingleInternalCore(inputData, singleModelName, false, out simData); - } - - /// - /// Simulate a single model to get the output including any additive inputs. - /// - /// - /// - /// - /// - public static bool SimulateSingle(TimeSeriesDataSet inputData, ISimulatableModel model, out TimeSeriesDataSet simData) - { - var plant = new PlantSimulator(new List { model }); - - return plant.Simulate(inputData, out simData ); - } + }*/ - /// - /// Simulates a single model for a unit dataset and adds the output to unitData.Y_meas of the unitData, optionally with noise - /// - /// the dataset to be simualted over, and where the Y_meas is updated with result - /// the model to be simulated - /// the amplitude of noise to be added to Y_meas - /// a seed value of the randm noise(specify so that tests are repeatable) - /// - public static bool SimulateSingleToYmeas(UnitDataSet unitData, ISimulatableModel model, double noiseAmplitude = 0, - int noiseSeed= 123) - { - (bool isOk, double[] y_proc) = SimulateSingleUnitDataWrapper(unitData, model); - - if (noiseAmplitude > 0) - { - // use a specific seed here, to avoid potential issues with "random unit tests" and not-repeatable - // errors. - Random rand = new Random(noiseSeed); - for (int k = 0; k < y_proc.Count(); k++) - { - y_proc[k] += (rand.NextDouble() - 0.5) * 2 * noiseAmplitude; - } - } - unitData.Y_meas = y_proc; - return isOk; - } - - /// - /// Simulates a single model for a unit dataset and adds the output to unitData.Y_meas of the unitData, optionally with noise - /// - /// the dataset to be simualted over, and where the Y_meas is updated with result - /// the model to be simulated - /// - public static (bool, double[]) SimulateSingleToYsim(UnitDataSet unitData, ISimulatableModel model) - { - (bool isOk, double[] y_proc) = SimulateSingleUnitDataWrapper(unitData, model); - unitData.Y_sim = y_proc; - return (isOk, y_proc); - } - - - /// - /// Simulates a single model given a unit data set - /// - /// - /// - /// - public static (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatableModel model) - { - return SimulateSingleUnitDataWrapper(unitData, model); - } - - /// - /// 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. - /// - /// contains a unit data set that must have U filled, Y_sim will be written here - /// model to simulate - /// a tuple, first aa true if able to simulate, otherwise false, second is the simulated time-series "y_proc" without any additive - static private (bool, double[]) SimulateSingleUnitDataWrapper(UnitDataSet unitData, ISimulatableModel model ) - { - const string defaultOutputName = "output"; - var inputData = new TimeSeriesDataSet(); - var singleModelName = "SimulateSingle"; - var modelCopy = model.Clone(singleModelName); - - if (unitData.Times != null) - inputData.SetTimeStamps(unitData.Times.ToList()); - else - { - inputData.CreateTimestamps(unitData.GetTimeBase()); - } - - var uNames = new List(); - for (int colIdx = 0; colIdx< unitData.U.GetNColumns(); colIdx++) - { - var uName = "U" + colIdx; - inputData.Add(uName, unitData.U.GetColumn(colIdx)); - uNames.Add(uName); - } - modelCopy.SetInputIDs(uNames.ToArray()); - { - inputData.Add(defaultOutputName, unitData.Y_meas); - modelCopy.SetOutputID(defaultOutputName); - } - - var sim = new PlantSimulator(new List { modelCopy }); - var isOk = sim.SimulateSingleInternalCore(inputData, singleModelName, false, out var simData); - - if(!isOk) - return (false, null); - - 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; - } - return (isOk, y_proc); - } - /// /// Simulate a single model(any ISimulatable model), using inputData as inputs, /// @@ -908,8 +693,18 @@ public bool Simulate (TimeSeriesDataSet inputData, out TimeSeriesDataSet simData var inputDataMinimal = new TimeSeriesDataSet(inputData); + // todo: disturbances could also instead be estimated in closed-loop? var didInit = init.ToSteadyStateAndEstimateDisturbances(ref inputDataMinimal, ref simData, compLoopDict); + // need to keep special track of pid-controlled outputs. + var pidControlledOutputsDict = DeterminePidControlledOutputs(); + // create internal "process outputs" for each such model + foreach (var pidOutput in pidControlledOutputsDict) + { + var signalID = pidOutput.Value; + simData.InitNewSignal(signalID, Double.NaN, N.Value); + } + if (!didInit) { Shared.GetParserObj().AddError("PlantSimulator failed to initalize."); @@ -965,6 +760,47 @@ public bool Simulate (TimeSeriesDataSet inputData, out TimeSeriesDataSet simData var model = modelDict[orderedSimulatorIDs.ElementAt(modelIdx)]; string[] inputIDs = model.GetBothKindsOfInputIDs(); + //////////////////////// + // before calculating a PID-model, treat the "junction" before it to get the right y_meas[k] = y_proc[k_1] + D[k] + if (model.GetProcessModelType() == ModelType.PID && timeIdx > 0) + { + var junctionSignalID = inputIDs.ElementAt((int)PidModelInputsIdx.Y_meas); + if (!pidControlledOutputsDict.ContainsKey(junctionSignalID)) + { + Shared.GetParserObj().AddError("PlantSimulator.Simulate() junction signal error \"" + junctionSignalID + "\""); + } + else + { + var modelID = pidControlledOutputsDict[junctionSignalID]; + if (modelID != null) // Simulataing a single PID-model not the whole loop + { + var value = simData.GetValue(modelID, timeIdx - 1); //y_proc[k-1] + if (modelID != null) + { + if (modelDict[modelID].GetAdditiveInputIDs() != null) // + D[k] + { + var additiveSignalsValues = GetValuesFromEitherDataset(modelDict[modelID], + modelDict[modelID].GetAdditiveInputIDs(), timeIdx, simData, inputDataMinimal); + foreach (var signalValue in additiveSignalsValues) + { + value += signalValue; + } + } + } + if (value.HasValue) + { + bool isOk = simData.AddDataPoint(junctionSignalID, timeIdx, value.Value); + } + else + { + Shared.GetParserObj().AddError("PlantSimulator.Simulate() error. Error calculating junction signal \"" + junctionSignalID + + "\""); + } + } + } + } + /////////////////////// + double[] inputVals = null; inputVals = GetValuesFromEitherDataset(model, inputIDs, lastGoodTimeIndex, simData, inputDataMinimal); if (inputVals == null) @@ -973,13 +809,30 @@ public bool Simulate (TimeSeriesDataSet inputData, out TimeSeriesDataSet simData "\" error retreiving input values."); return false; } + double[] outputVal = model.Iterate(inputVals, timeBase_s); - bool isOk = simData.AddDataPoint(model.GetOutputID(),timeIdx,outputVal[0]); - if (!isOk) + + if (pidControlledOutputsDict.Keys.Contains(model.GetOutputID())) { - Shared.GetParserObj().AddError("PlantSimulator.Simulate() failed. Unable to add data point for \"" - + model.GetOutputID() + "\", indicating an error in initalization. "); - return false; + // for pid-controlled outputs, save the result with the name of the process ID, to be used in the + // next iteration, possibly in combination with a disturbance signal + bool isOk = simData.AddDataPoint(model.GetID(), timeIdx, outputVal[0]); + if (!isOk) + { + Shared.GetParserObj().AddError("PlantSimulator.Simulate() failed. Unable to add data point for \"" + + model.GetOutputID() + "\", indicating an error in initalization. "); + return false; + } + } + else + { + bool isOk = simData.AddDataPoint(model.GetOutputID(), timeIdx, outputVal[0]); + if (!isOk) + { + Shared.GetParserObj().AddError("PlantSimulator.Simulate() failed. Unable to add data point for \"" + + model.GetOutputID() + "\", indicating an error in initalization. "); + return false; + } } if (outputVal.Length > 1) { @@ -993,6 +846,8 @@ public bool Simulate (TimeSeriesDataSet inputData, out TimeSeriesDataSet simData } } } + if (inputDataMinimal != null) + if (inputDataMinimal.GetTimeStamps() != null) simData.SetTimeStamps(inputDataMinimal.GetTimeStamps().ToList()); PlantFitScore = FitScoreCalculator.GetPlantWideSimulated(this, inputData, simData); diff --git a/Dynamic/PlantSimulator/PlantSimulatorHelper.cs b/Dynamic/PlantSimulator/PlantSimulatorHelper.cs new file mode 100644 index 00000000..d488eb60 --- /dev/null +++ b/Dynamic/PlantSimulator/PlantSimulatorHelper.cs @@ -0,0 +1,256 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace TimeSeriesAnalysis.Dynamic +{ + /// + /// Convenience functions for using PlantSimulator + /// + public class PlantSimulatorHelper + { + /// + /// Create a PlantSimulator and TimeSeriesDataSet from a UnitDataSet, PidModel and UnitModel to do closed-loop simulations + /// + /// The feedback loop has NO disturbance signal added, but this can be added to the returned PlantSimulator as needed. + /// + /// + /// + /// + /// + /// + /// a simulator object and a dataset object that is ready to be simulated with Simulate() + public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoopNoDisturbance(UnitDataSet unitDataSet, PidModel pidModel, + UnitModel unitModel, int pidInputIdx = 0) + { + var plantSim = new PlantSimulator( + new List { pidModel, unitModel }); + var signalId1 = plantSim.ConnectModels(unitModel, pidModel); + var signalId2 = plantSim.ConnectModels(pidModel, unitModel, pidInputIdx); + + var inputData = new TimeSeriesDataSet(); + inputData.Add(signalId1, (double[])unitDataSet.Y_meas.Clone()); + + for (int curColIdx = 0; curColIdx < unitDataSet.U.GetNColumns(); curColIdx++) + { + if (curColIdx == pidInputIdx) + { + inputData.Add(signalId2, (double[])unitDataSet.U.GetColumn(pidInputIdx).Clone()); + } + else + { + inputData.Add(plantSim.AddExternalSignal(unitModel, SignalType.External_U, curColIdx), + (double[])unitDataSet.U.GetColumn(curColIdx).Clone()); + } + } + inputData.Add(plantSim.AddExternalSignal(pidModel, SignalType.Setpoint_Yset), (double[])unitDataSet.Y_setpoint.Clone()); + inputData.CreateTimestamps(unitDataSet.GetTimeBase()); + inputData.SetIndicesToIgnore(unitDataSet.IndicesToIgnore); + return (plantSim, inputData); + } + + /// + /// Create a feedback loop, where the process model has an additive disturbance that is to be estimated. + /// + /// + /// + /// + /// + /// a simulator object and a dataset object that is ready to be simulated with Simulate() + public static (PlantSimulator, TimeSeriesDataSet) CreateFeedbackLoopWithEstimatedDisturbance(UnitDataSet unitDataSet, PidModel pidModel, + UnitModel unitModel, int pidInputIdx = 0) + { + // vital that signal follows naming convention, otherwise it will not be estimated, but should be provided. + unitModel.AddSignalToOutput(SignalNamer.EstDisturbance(unitModel)); + (var sim, var data) = CreateFeedbackLoopNoDisturbance(unitDataSet, pidModel, unitModel, pidInputIdx); + return (sim, data); + } + + + /// + /// Returns a unit data set for a given UnitModel. + /// + /// + /// + /// + /// a tuple with a bool indicating if it was a success as item1, and the dataset as item2 + public static (bool,UnitDataSet) GetUnitDataSetForLoop(TimeSeriesDataSet inputData, PidModel pidModel, UnitModel unitModel) + { + UnitDataSet dataset = new UnitDataSet(); + var inputIDs = unitModel.GetModelInputIDs(); + dataset.U = new double[inputData.GetLength().Value, inputIDs.Length]; + bool success = true; + dataset.Times = inputData.GetTimeStamps(); + + var outputID = unitModel.GetOutputID(); + if (inputData.ContainsSignal(outputID)) + dataset.Y_meas = inputData.GetValues(outputID); + else + success = false; + for (int inputIDidx = 0; inputIDidx < inputIDs.Length; inputIDidx++) + { + var inputID = inputIDs[inputIDidx]; + if (inputData.ContainsSignal(inputID)) + { + var curCol = inputData.GetValues(inputID); + dataset.U.WriteColumn(inputIDidx, curCol); + } + else + { + success = false; + } + } + + var setpointID = pidModel.GetModelInputIDs().ElementAt((int)PidModelInputsIdx.Y_setpoint); + if (inputData.ContainsSignal(setpointID)) + dataset.Y_setpoint = inputData.GetValues(setpointID); + else + success = false; + + return (success,dataset); + } + + /// + /// Simulate a single model to get the output including any additive inputs. + /// + /// + /// + /// + /// + public static bool SimulateSingle(TimeSeriesDataSet inputData, ISimulatableModel model, out TimeSeriesDataSet simData) + { + PlantSimulator plant = null; + /* if (model.GetProcessModelType() == ModelType.SubProcess) + { + var modelClone = (UnitModel)model.Clone("clone"); + modelClone.RemoveAdditiveInputs(); + + plant = new PlantSimulator(new List { modelClone }); + } + else*/ + { + plant = new PlantSimulator(new List { model }); + } + return plant.Simulate(inputData, out simData); + } + + + /// + /// Simulates a single model given a unit data set + /// + /// + /// + /// + public static (bool, double[]) SimulateSingle(UnitDataSet unitData, ISimulatableModel model) + { + return SimulateSingleUnitDataWrapper(unitData, model); + } + + /// + /// Simulates a single model for a unit dataset and adds the output to unitData.Y_meas of the unitData, optionally with noise + /// + /// the dataset to be simualted over, and where the Y_meas is updated with result + /// the model to be simulated + /// the amplitude of noise to be added to Y_meas + /// a seed value of the randm noise(specify so that tests are repeatable) + /// + public static bool SimulateSingleToYmeas(UnitDataSet unitData, ISimulatableModel model, double noiseAmplitude = 0, + int noiseSeed = 123) + { + (bool isOk, double[] y_proc) = SimulateSingleUnitDataWrapper(unitData, model); + + if (noiseAmplitude > 0) + { + // use a specific seed here, to avoid potential issues with "random unit tests" and not-repeatable + // errors. + Random rand = new Random(noiseSeed); + for (int k = 0; k < y_proc.Count(); k++) + { + y_proc[k] += (rand.NextDouble() - 0.5) * 2 * noiseAmplitude; + } + } + unitData.Y_meas = y_proc; + return isOk; + } + + /// + /// Simulates a single model for a unit dataset and adds the output to unitData.Y_meas of the unitData, optionally with noise + /// + /// the dataset to be simualted over, and where the Y_meas is updated with result + /// the model to be simulated + /// + public static (bool, double[]) SimulateSingleToYsim(UnitDataSet unitData, ISimulatableModel model) + { + (bool isOk, double[] y_proc) = SimulateSingleUnitDataWrapper(unitData, model); + unitData.Y_sim = y_proc; + return (isOk, y_proc); + } + + + /// + /// 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. + /// + /// contains a unit data set that must have U filled, Y_sim will be written here + /// model to simulate + /// a tuple, first aa true if able to simulate, otherwise false, second is the simulated time-series "y_proc" without any additive + static private (bool, double[]) SimulateSingleUnitDataWrapper(UnitDataSet unitData, ISimulatableModel model) + { + const string defaultOutputName = "output"; + var inputData = new TimeSeriesDataSet(); + var singleModelName = "SimulateSingle"; + var modelCopy = model.Clone(singleModelName); + + if (unitData.Times != null) + inputData.SetTimeStamps(unitData.Times.ToList()); + else + { + inputData.CreateTimestamps(unitData.GetTimeBase()); + } + + var uNames = new List(); + for (int colIdx = 0; colIdx < unitData.U.GetNColumns(); colIdx++) + { + var uName = "U" + colIdx; + inputData.Add(uName, unitData.U.GetColumn(colIdx)); + uNames.Add(uName); + } + modelCopy.SetInputIDs(uNames.ToArray()); + { + inputData.Add(defaultOutputName, unitData.Y_meas); + modelCopy.SetOutputID(defaultOutputName); + } + + var isOk = PlantSimulatorHelper.SimulateSingle(inputData, modelCopy, out var simData); + + // var sim = new PlantSimulator(new List { modelCopy }); + // var isOk = sim.Simulate(inputData, out var simData); + + if (!isOk) + return (false, null); + + 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; + } + return (isOk, y_proc); + } + + + + } +} diff --git a/Dynamic/PlantSimulator/PlantSimulatorInitalizer.cs b/Dynamic/PlantSimulator/PlantSimulatorInitalizer.cs index 612a430a..027f5d70 100644 --- a/Dynamic/PlantSimulator/PlantSimulatorInitalizer.cs +++ b/Dynamic/PlantSimulator/PlantSimulatorInitalizer.cs @@ -114,7 +114,7 @@ public bool ToSteadyStateAndEstimateDisturbances(ref TimeSeriesDataSet inputData { simulatedSignals.Add(simulator.modelDict[modelName].GetOutputID()); } - // in addtion, any signal addted to signalValuesAtT0, but not in inputData is to be simulated. + // in addtion, any signal added to signalValuesAtT0, but not in inputData is to be simulated. // this includes disturbances foreach (string signalID in signalValuesAtT0.Keys) { @@ -240,13 +240,14 @@ private bool InitComputationalLoops(Dictionary> compLoopDic /// For a plant, go through and find each plant/pid-controller and attempt to estimate the disturbance. /// For the disturbance to be estimateable,the inputs "u_meas" and the outputs "y_meas" for each "process" in /// each pid-process loop needs to be given in inputData. - /// The estimated disturbance signal is addes to simData + /// The estimated disturbance signal is added to simData /// /// note that for closed loop systems u and y signals are removed(these are used to estimate disturbance, removing them triggers code in PlantSimulator to re-estimate them) /// /// /// true if everything went ok, otherwise false - private bool EstimateDisturbances(ref TimeSeriesDataSet inputData, ref TimeSeriesDataSet simData,ref Dictionary signalValuesAtT0) + private bool EstimateDisturbances(ref TimeSeriesDataSet inputData, ref TimeSeriesDataSet simData, + ref Dictionary signalValuesAtT0) { // find all PID-controllers List pidIDs = new List(); @@ -257,6 +258,7 @@ private bool EstimateDisturbances(ref TimeSeriesDataSet inputData, ref TimeSerie pidIDs.Add(model.Key); } } + foreach (var pidID in pidIDs) { var upstreamModels = simulator.connections.GetAllUpstreamModels(pidID); @@ -264,49 +266,79 @@ private bool EstimateDisturbances(ref TimeSeriesDataSet inputData, ref TimeSerie if (upstreamModels.Count == 0) continue; var processId = upstreamModels.First(); - - var isOK = simulator.SimulateSingleWithoutAdditive(inputData, processId, - out TimeSeriesDataSet singleSimDataSetWithDisturbance); - if (isOK) + var estDisturbanceId = SignalNamer.EstDisturbance(processId); + + if (inputData.ContainsSignal(estDisturbanceId)) + continue; + + bool doV1 = true; + + if (doV1) { - var estDisturbanceId = SignalNamer.EstDisturbance(processId); - if (singleSimDataSetWithDisturbance.ContainsSignal(estDisturbanceId)) + + var isOK = simulator.SimulateSingleWithoutAdditive(inputData, processId, + out var singleSimDataSetWithDisturbance); + if (isOK) { - var estDisturbance = singleSimDataSetWithDisturbance.GetValues(estDisturbanceId); - if (estDisturbance == null) - continue; - if ((new Vec()).IsAllNaN(estDisturbance)) - continue; - // add signal if everything is ok. - simData.Add(estDisturbanceId, estDisturbance); - // todo: remove pid input pid-u and output y from inputdata(we want to re-estimate it, we have used it to estimate the disturbance) - // an alterntive to this would hav been to to alter the code in the plant simulator to add signals in simData output that are duplicates of signal names in inputData - // or better: make an internal "stripped" version of "inputData" + + if (singleSimDataSetWithDisturbance.ContainsSignal(estDisturbanceId)) { - { - string y_meas_signal = simulator.modelDict[processId].GetOutputID(); - double? value = inputData.GetValue(y_meas_signal, 0); - if (value.HasValue) - { - signalValuesAtT0.Add(y_meas_signal, value.Value); - inputData.Remove(y_meas_signal); - } - else - return false; - } + var estDisturbance = singleSimDataSetWithDisturbance.GetValues(estDisturbanceId); + if (estDisturbance == null) + continue; + if ((new Vec()).IsAllNaN(estDisturbance)) + continue; + // add signal if everything is ok. + simData.Add(estDisturbanceId, estDisturbance); + } + + } + } + else //: new version that uses DisturbanceCalculator + { + var processModel = (UnitModel)simulator.modelDict[processId]; + var pidModel = (PidModel)simulator.modelDict[pidID]; + var pidParams = pidModel.GetModelParameters(); + var pidOutName = pidModel.outputID; + var pidInputIdx = 0; + (var isDataOk,var dataset) = PlantSimulatorHelper.GetUnitDataSetForLoop(inputData, pidModel, processModel); + if (isDataOk) + { + if (dataset.U.Length > 1) + { + string pidOutId = pidModel.outputID; + for (int i = 0; i < processModel.ModelInputIDs.Length; i++) { - string u_pid_signal = simulator.modelDict[pidID].GetOutputID(); - double? value = inputData.GetValue(u_pid_signal, 0); - if (value.HasValue) - { - signalValuesAtT0.Add(u_pid_signal,value.Value); - inputData.Remove(u_pid_signal); - } - else - return false; + if (processModel.ModelInputIDs[i] == pidOutId) + pidInputIdx = i; } } + var distResult = DisturbanceCalculator.CalculateDisturbanceVector(dataset, processModel, pidInputIdx, pidParams); + simData.Add(estDisturbanceId, distResult.d_est); + } + } + + { + string y_meas_signal = simulator.modelDict[processId].GetOutputID(); + double? value = inputData.GetValue(y_meas_signal, 0); + if (value.HasValue) + { + signalValuesAtT0.Add(y_meas_signal, value.Value); + inputData.Remove(y_meas_signal); + } + else + return false; + } + { + string u_pid_signal = simulator.modelDict[pidID].GetOutputID(); + double? value = inputData.GetValue(u_pid_signal, 0); + if (value.HasValue) + { + signalValuesAtT0.Add(u_pid_signal, value.Value); + inputData.Remove(u_pid_signal); } + else + return false; } } return true; diff --git a/Dynamic/SimulatableModels/UnitModel.cs b/Dynamic/SimulatableModels/UnitModel.cs index 5d3fef75..25f076a2 100644 --- a/Dynamic/SimulatableModels/UnitModel.cs +++ b/Dynamic/SimulatableModels/UnitModel.cs @@ -42,12 +42,12 @@ namespace TimeSeriesAnalysis.Dynamic /// /// See also: /// - public class UnitModel : ModelBaseClass, ISimulatableModel + public class UnitModel : ModelBaseClass, ISimulatableModel { /// /// The parameters of the UnitModel. /// - public UnitParameters modelParameters; + public UnitParameters modelParameters; // first order dynamics : private LowPass lowPass1order; @@ -56,21 +56,21 @@ public class UnitModel : ModelBaseClass, ISimulatableModel private SecondOrder secondOrder; - + // time-delay private TimeDelay delayObj; private bool isFirstIteration; private double[] lastGoodValuesOfInputs; - private UnitDataSet FittedDataSet=null; + private UnitDataSet FittedDataSet = null; private List TimeDelayEstWarnings { get; } /// /// Empty constructor /// - public UnitModel(){} + public UnitModel() { } /// /// Constructor @@ -78,7 +78,7 @@ public UnitModel(){} /// model parameter object /// a unique string that identifies this model in larger process models [JsonConstructor] - public UnitModel(UnitParameters modelParameters, string ID="not_named") + public UnitModel(UnitParameters modelParameters, string ID = "not_named") { processModelType = ModelType.SubProcess; this.ID = ID; @@ -90,7 +90,7 @@ public UnitModel(UnitParameters modelParameters, string ID="not_named") /// /// /// a unique string that identifies this model in larger process models - public UnitModel(UnitParameters modelParameters, UnitDataSet dataSet,string ID = "not_named") + public UnitModel(UnitParameters modelParameters, UnitDataSet dataSet, string ID = "not_named") { processModelType = ModelType.SubProcess; this.ID = ID; @@ -110,16 +110,16 @@ public bool IsModelSimulatable(out string explainStr) explainStr = "modelParameters is null"; return false; } - /* if (modelParameters.LinearGains == null) - { - explainStr = "LinearGains is null"; - return false; - } - if (modelParameters.LinearGains.Length == 0) - { - explainStr = "LinearGains is empty"; - return false; - }*/ + /* if (modelParameters.LinearGains == null) + { + explainStr = "LinearGains is null"; + return false; + } + if (modelParameters.LinearGains.Length == 0) + { + explainStr = "LinearGains is empty"; + return false; + }*/ if (ModelInputIDs == null) { explainStr = "ModelInputIDs is null"; @@ -162,9 +162,9 @@ public void InitSim(UnitParameters modelParameters) } if (doSim) { - this.lastGoodValuesOfInputs = Vec.Fill(Double.NaN, + this.lastGoodValuesOfInputs = Vec.Fill(Double.NaN, GetLengthOfInputVector()); - this.SetInputIDs(new string[GetLengthOfInputVector()]); + this.SetInputIDs(new string[GetLengthOfInputVector()]); } } @@ -212,7 +212,7 @@ override public int GetLengthOfInputVector() return inputIDs.Length; } } - + /// /// Get the object of model parameters contained in the model. /// @@ -241,9 +241,9 @@ public void SetModelParameters(UnitParameters parameters) /// /// /// - public double? GetSteadyStateInput(double x0, int inputIdx=0, double[] givenInputs=null) + public double? GetSteadyStateInput(double x0, int inputIdx = 0, double[] givenInputs = null) { - double u0=0; + double u0 = 0; if (modelParameters.LinearGains == null) { @@ -300,7 +300,7 @@ public void SetModelParameters(UnitParameters parameters) u0 = 0; if (modelParameters.U0 != null) { - u0 += modelParameters.U0[inputIdx]; + u0 += modelParameters.U0[inputIdx]; } u0 += y_contributionFromInput / modelParameters.LinearGains[inputIdx]; @@ -315,18 +315,18 @@ public void SetModelParameters(UnitParameters parameters) double b = modelParameters.LinearGains[inputIdx]; double c = -y_contributionFromInput; double[] quadSolution = SolveQuadratic(a, b, c); - double chosenU=0; + double chosenU = 0; if (quadSolution.Length == 1) - { - chosenU = quadSolution[0]; + { + chosenU = quadSolution[0]; } else if (quadSolution.Length == 2) - { + { chosenU = Math.Min(quadSolution[0], quadSolution[1]); } if (modelParameters.U0 == null) { - u0 = chosenU ; + u0 = chosenU; } else { @@ -348,7 +348,7 @@ private double CalculateLinearProcessGainTerm(int inputIndex, double u) double processGainTerm = 0; if (modelParameters.U0 != null) { - processGainTerm += modelParameters.LinearGains[inputIndex] * (u- modelParameters.U0[inputIndex]); + processGainTerm += modelParameters.LinearGains[inputIndex] * (u - modelParameters.U0[inputIndex]); } else { @@ -388,12 +388,12 @@ private double CalculateCurvatureProcessGainTerm(int inputIndex, double u) if (modelParameters.U0 != null) { curvatureTerm += modelParameters.Curvatures[inputIndex] * - Math.Pow((u - modelParameters.U0[inputIndex]) , 2) / uNorm; + Math.Pow((u - modelParameters.U0[inputIndex]), 2) / uNorm; } else { curvatureTerm += modelParameters.Curvatures[inputIndex] * - Math.Pow(u , 2) / uNorm; + Math.Pow(u, 2) / uNorm; } } return curvatureTerm; @@ -405,7 +405,7 @@ private double CalculateCurvatureProcessGainTerm(int inputIndex, double u) /// /// /// - private double CalculateSteadyStateWithoutAdditive(double[] inputs, double badValueIndicator=-9999) + private double CalculateSteadyStateWithoutAdditive(double[] inputs, double badValueIndicator = -9999) { double x_ss = modelParameters.Bias; // inputs U may include a disturbance as the last entry @@ -439,7 +439,7 @@ private double CalculateSteadyStateWithoutAdditive(double[] inputs, double badVa /// vector of input values /// /// - public double? GetSteadyStateOutput(double[] u, double badDataID=-9999) + public double? GetSteadyStateOutput(double[] u, double badDataID = -9999) { if (modelParameters.LinearGains == null) return 0; @@ -478,9 +478,9 @@ override public SignalType GetOutputSignalType() /// NaN is returned if model was not able to be identfied, or if no good U values have been given yet. /// If some data points in U inputs are NaN or equal to badValueIndicator, the last good value is returned /// - public double[] Iterate(double[] inputs, double timeBase_s,double badValueIndicator=-9999) + public double[] Iterate(double[] inputs, double timeBase_s, double badValueIndicator = -9999) { - if (modelParameters.Fitting!= null) + if (modelParameters.Fitting != null) { if (!modelParameters.Fitting.WasAbleToIdentify) { @@ -500,14 +500,14 @@ public double[] Iterate(double[] inputs, double timeBase_s,double badValueIndica // this is the case for newly created models that have not yet been fitted if (modelParameters.LinearGains == null) { - return new double[] { 0 }; + return new double[] { 0 }; } // notice! the model does not use the parameters [a,b,c,unorm,td.u0] to simulate the process model // instead it calculates the steady-state and then filters the steady-state with LowPass to get the appropriate time constant // - so it uses the parameters [linearGain, curvatureGain, Timeconstant,td] - double x_ss = CalculateSteadyStateWithoutAdditive(inputs,badValueIndicator); + double x_ss = CalculateSteadyStateWithoutAdditive(inputs, badValueIndicator); // nb! if first iteration, start model at steady-state double x_dynamic = x_ss; @@ -528,7 +528,7 @@ public double[] Iterate(double[] inputs, double timeBase_s,double badValueIndica double y = 0; if (modelParameters.TimeDelay_s <= 0) { - y = x_dynamic; + y = x_dynamic; } else { @@ -556,10 +556,10 @@ public double[] Iterate(double[] inputs, double timeBase_s,double badValueIndica } } if (y_internal.HasValue) - return new double[] { y, y_internal.Value}; + return new double[] { y, y_internal.Value }; else return new double[] { y }; - } + } /// /// Is the model static or dynamic? @@ -567,7 +567,7 @@ public double[] Iterate(double[] inputs, double timeBase_s,double badValueIndica /// Returns true if the model is static(no time constant or time delay terms),otherwise false. public bool IsModelStatic() { - return modelParameters.TimeConstant_s == 0 && modelParameters.TimeDelay_s == 0; + return modelParameters.TimeConstant_s == 0 && modelParameters.TimeDelay_s == 0; } @@ -604,6 +604,14 @@ private static double[] SolveQuadratic(double a, double b, double c) } } + /// + /// Removes the addtive inputs. + /// + public void RemoveAdditiveInputs() + { + additiveInputIDs = new List(); + } + /// /// Create a nice human-readable summary of all the important data contained in the model object. diff --git a/TimeSeriesAnalysis.Tests/Examples/ProcessControl.cs b/TimeSeriesAnalysis.Tests/Examples/ProcessControl.cs index 4bf5505d..20194458 100644 --- a/TimeSeriesAnalysis.Tests/Examples/ProcessControl.cs +++ b/TimeSeriesAnalysis.Tests/Examples/ProcessControl.cs @@ -160,7 +160,7 @@ var disturbanceModel inputData.Add(simNoFeedF.AddExternalSignal(pidModel, SignalType.Setpoint_Yset), TimeSeriesCreator.Constant(60, 600)); inputData.Add(simNoFeedF.AddExternalSignal(disturbanceModel, SignalType.External_U), - TimeSeriesCreator.Step(300, 600, 25, 0)); + TimeSeriesCreator.Step(100, 600, 25, 0)); inputData.CreateTimestamps(timeBase_s); var isOk = simNoFeedF.Simulate(inputData,out var dataNoFeedF); @@ -172,7 +172,7 @@ var disturbanceModel dataNoFeedF.GetValues(pidModel.GetID(),SignalType.PID_U), inputData.GetValues(disturbanceModel.GetID(),SignalType.External_U) }, - new List { "y1=y_run1", "y1=y_setpoint", "y2=y_dist[right]", "y3=u_pid", "y3=u_dist" }, timeBase_s, "FeedForwardEx1"); + new List { "y1=y_meas", "y1=y_setpoint", "y2=y_dist[right]", "y3=u_pid", "y3=u_dist" }, timeBase_s, "FeedForwardEx1"); #endregion @@ -527,7 +527,6 @@ public void IntegralOscillations(double KP) // unless the disturbance oscillates // integral oscillations are often thought to be created by - double noiseAmp = 1; UnitParameters modelParameters = new UnitParameters { diff --git a/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs b/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs index c621e55e..86404f49 100644 --- a/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs +++ b/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs @@ -68,7 +68,7 @@ public void NonlinearUnitModel() var refData = new UnitDataSet(); refData.U = U; refData.CreateTimeStamps(timeBase_s); - PlantSimulator.SimulateSingleToYsim(refData, refModel); + PlantSimulatorHelper.SimulateSingleToYsim(refData, refModel); // simulate the nonlinear model UnitParameters designParameters = new UnitParameters @@ -85,7 +85,7 @@ public void NonlinearUnitModel() var idDataSet = new UnitDataSet(); idDataSet.U = U; idDataSet.CreateTimeStamps(timeBase_s); - PlantSimulator.SimulateSingleToYmeas(idDataSet, trueModel,noiseAmplitude); + PlantSimulatorHelper.SimulateSingleToYmeas(idDataSet, trueModel,noiseAmplitude); // do identification of unit model FittingSpecs fittingSpecs = new FittingSpecs(designParameters.U0, designParameters.UNorm); diff --git a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_DynamicSISO.cs b/TimeSeriesAnalysis.Tests/Test/DisturbanceID/ClosedLoopIdentifierTester_DynamicSISO.cs similarity index 100% rename from TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_DynamicSISO.cs rename to TimeSeriesAnalysis.Tests/Test/DisturbanceID/ClosedLoopIdentifierTester_DynamicSISO.cs diff --git a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_MISO.cs b/TimeSeriesAnalysis.Tests/Test/DisturbanceID/ClosedLoopIdentifierTester_MISO.cs similarity index 93% rename from TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_MISO.cs rename to TimeSeriesAnalysis.Tests/Test/DisturbanceID/ClosedLoopIdentifierTester_MISO.cs index 4844c807..eb774fdf 100644 --- a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_MISO.cs +++ b/TimeSeriesAnalysis.Tests/Test/DisturbanceID/ClosedLoopIdentifierTester_MISO.cs @@ -26,6 +26,8 @@ class ClosedLoopIdentifierTester_MISO Bias = 10 }; + + UnitParameters dynamicModelParameters = new UnitParameters { TimeConstant_s = 8, @@ -61,7 +63,7 @@ public void CommonPlotAndAsserts(UnitDataSet pidDataSet, double[] estDisturbance Assert.IsTrue(estDisturbance != null); string caseId = TestContext.CurrentContext.Test.Name.Replace("(", "_"). Replace(")", "_").Replace(",", "_") + "y"; - bool doDebugPlot = false; + bool doDebugPlot = true; if (doDebugPlot) { Shared.EnablePlots(); @@ -99,27 +101,27 @@ public void CommonPlotAndAsserts(UnitDataSet pidDataSet, double[] estDisturbance } } - [TestCase(0,false)] - [TestCase(1,false)] + + [TestCase(0, false)] [TestCase(0, true)] - public void StaticMISO_SetpointChanges_no_disturbance_detectsProcessOk(int pidInputIdx, bool doNegative) + public void Static2Input_SetpointChanges_NOdisturbance_detectsProcessOk(int pidInputIdx, bool doNegative) { var trueDisturbance = TimeSeriesCreator.Constant(0, N); - var externalU1 = TimeSeriesCreator.Step(N/8, N, 5, 10); - var externalU2 = TimeSeriesCreator.Step(N*5/8, N, 2, 1); - var yset = TimeSeriesCreator.Step(N*3/8, N, 20, 18); - var trueParamters = staticModelParameters.CreateCopy(); - if (pidInputIdx == 1) + var externalU1 = TimeSeriesCreator.Step(N / 8, N, 15, 20); + var yset = TimeSeriesCreator.Step(N * 3 / 8, N, 20, 18); + UnitParameters twoInputModelParameters = new UnitParameters { - double[] oldGains = new double[3]; - trueParamters.LinearGains.CopyTo(oldGains); - trueParamters.LinearGains = new double[] { oldGains[1], oldGains[0], oldGains[2] }; - } - GenericMISODisturbanceTest(new UnitModel(trueParamters, "StaticProcess"), trueDisturbance, externalU1,externalU2, - doNegative, true,yset, pidInputIdx,10,false, isStatic); + TimeConstant_s = 0, + LinearGains = new double[] { 0.5, 0.25 }, + TimeDelay_s = 0, + Bias = 10 + }; + GenericMISODisturbanceTest(new UnitModel(twoInputModelParameters, "StaticTwoInputsProcess"), + trueDisturbance, externalU1, null, doNegative, true, yset, pidInputIdx, 10, false, isStatic); } + [TestCase(0, false, false), NonParallelizable] [TestCase(1, false, false)] // when a third input is added, estimation seems to fail. @@ -229,7 +231,7 @@ public void DynamicMISO_externalUchanges_NoDisturbance_NOsetpointchange_DetectsP { TimeConstant_s = 10, LinearGains = new double[] { 1, 0.5}, - TimeDelay_s = 5, + TimeDelay_s = 0, Bias = 5 }; @@ -325,6 +327,7 @@ public void GenericMISODisturbanceTest (UnitModel trueProcessModel, double[] tr pidDataSet.Y_setpoint[50] = Double.NaN; pidDataSet.Y_meas[400] = Double.NaN; pidDataSet.U[500,0] = Double.NaN; + Console.WriteLine("BAD DATA ADDED!!!"); } (var identifiedModel, var estDisturbance) = ClosedLoopUnitIdentifier.Identify(pidDataSet, pidModel1.GetModelParameters(), pidInputIdx); diff --git a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs b/TimeSeriesAnalysis.Tests/Test/DisturbanceID/ClosedLoopIdentifierTester_StaticSISO.cs similarity index 100% rename from TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs rename to TimeSeriesAnalysis.Tests/Test/DisturbanceID/ClosedLoopIdentifierTester_StaticSISO.cs diff --git a/TimeSeriesAnalysis.Tests/Tests/CluiCommonTests.cs b/TimeSeriesAnalysis.Tests/Test/DisturbanceID/CluiCommonTests.cs similarity index 100% rename from TimeSeriesAnalysis.Tests/Tests/CluiCommonTests.cs rename to TimeSeriesAnalysis.Tests/Test/DisturbanceID/CluiCommonTests.cs diff --git a/TimeSeriesAnalysis.Tests/Tests/DisturbanceCalculatorTests.cs b/TimeSeriesAnalysis.Tests/Test/DisturbanceID/DisturbanceCalculatorTests.cs similarity index 100% rename from TimeSeriesAnalysis.Tests/Tests/DisturbanceCalculatorTests.cs rename to TimeSeriesAnalysis.Tests/Test/DisturbanceID/DisturbanceCalculatorTests.cs diff --git a/TimeSeriesAnalysis.Tests/Tests/Array2DUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/Fundamentals/Array2DUnitTests.cs similarity index 98% rename from TimeSeriesAnalysis.Tests/Tests/Array2DUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/Fundamentals/Array2DUnitTests.cs index 36bd1858..5e053d84 100644 --- a/TimeSeriesAnalysis.Tests/Tests/Array2DUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/Fundamentals/Array2DUnitTests.cs @@ -6,7 +6,7 @@ using NUnit.Framework; using TimeSeriesAnalysis; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Fundamentals { [TestFixture] class ArrayUnitTests diff --git a/TimeSeriesAnalysis.Tests/Tests/IndexTests.cs b/TimeSeriesAnalysis.Tests/Test/Fundamentals/IndexTests.cs similarity index 96% rename from TimeSeriesAnalysis.Tests/Tests/IndexTests.cs rename to TimeSeriesAnalysis.Tests/Test/Fundamentals/IndexTests.cs index 8f0d32a3..9724114c 100644 --- a/TimeSeriesAnalysis.Tests/Tests/IndexTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/Fundamentals/IndexTests.cs @@ -7,7 +7,7 @@ using TimeSeriesAnalysis.Utility; using NUnit.Framework; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Fundamentals { [TestFixture] class IndexTest diff --git a/TimeSeriesAnalysis.Tests/Tests/LowPassUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/Fundamentals/LowPassUnitTests.cs similarity index 97% rename from TimeSeriesAnalysis.Tests/Tests/LowPassUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/Fundamentals/LowPassUnitTests.cs index dc530e4b..6b7530ab 100644 --- a/TimeSeriesAnalysis.Tests/Tests/LowPassUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/Fundamentals/LowPassUnitTests.cs @@ -8,7 +8,7 @@ using TimeSeriesAnalysis.Dynamic; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Fundamentals { [TestFixture] class FilterUnitTests diff --git a/TimeSeriesAnalysis.Tests/Tests/MatrixUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/Fundamentals/MatrixUnitTests.cs similarity index 98% rename from TimeSeriesAnalysis.Tests/Tests/MatrixUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/Fundamentals/MatrixUnitTests.cs index c01dd67c..4017e3e7 100644 --- a/TimeSeriesAnalysis.Tests/Tests/MatrixUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/Fundamentals/MatrixUnitTests.cs @@ -6,7 +6,7 @@ using NUnit.Framework; using TimeSeriesAnalysis; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Fundamentals { [TestFixture] class MatrixUnitTests diff --git a/TimeSeriesAnalysis.Tests/Tests/SignificantDigitsTests.cs b/TimeSeriesAnalysis.Tests/Test/Fundamentals/SignificantDigitsTests.cs similarity index 96% rename from TimeSeriesAnalysis.Tests/Tests/SignificantDigitsTests.cs rename to TimeSeriesAnalysis.Tests/Test/Fundamentals/SignificantDigitsTests.cs index 46f97111..44a4937d 100644 --- a/TimeSeriesAnalysis.Tests/Tests/SignificantDigitsTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/Fundamentals/SignificantDigitsTests.cs @@ -8,7 +8,7 @@ using TimeSeriesAnalysis.Utility; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Fundamentals { [TestFixture] internal class SignificantDigitsTests diff --git a/TimeSeriesAnalysis.Tests/Tests/VecExtensionsMethodsUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/Fundamentals/VecExtensionsMethodsUnitTests.cs similarity index 89% rename from TimeSeriesAnalysis.Tests/Tests/VecExtensionsMethodsUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/Fundamentals/VecExtensionsMethodsUnitTests.cs index 4f3a16b7..ecb071da 100644 --- a/TimeSeriesAnalysis.Tests/Tests/VecExtensionsMethodsUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/Fundamentals/VecExtensionsMethodsUnitTests.cs @@ -6,7 +6,7 @@ using NUnit.Framework; using TimeSeriesAnalysis; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Fundamentals { [TestFixture] class VecExtensionsUnitTests diff --git a/TimeSeriesAnalysis.Tests/Tests/VecUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/Fundamentals/VecTests.cs similarity index 99% rename from TimeSeriesAnalysis.Tests/Tests/VecUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/Fundamentals/VecTests.cs index d5ba65ef..20e7acd9 100644 --- a/TimeSeriesAnalysis.Tests/Tests/VecUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/Fundamentals/VecTests.cs @@ -5,10 +5,10 @@ using TimeSeriesAnalysis.Utility; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Fundamentals { [TestFixture] - class VecUnitTests + class VecTests { [TestCase(new int[] { 1, 2, 3 }, new int[] { 1, 2, 3, 4 })] diff --git a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorMISOTests.cs b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/AdvancedControlExamples.cs similarity index 84% rename from TimeSeriesAnalysis.Tests/Tests/PlantSimulatorMISOTests.cs rename to TimeSeriesAnalysis.Tests/Test/PlantSimulations/AdvancedControlExamples.cs index c880bdcd..1f4947ef 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorMISOTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/AdvancedControlExamples.cs @@ -20,28 +20,31 @@ namespace TimeSeriesAnalysis.Test.PlantSimulations /// /// [TestFixture] - class Control + class AdvancedControlExamples { [Test] public void CascadeControl() { - Shared.DisablePlots(); - + // Shared.DisablePlots(); + // Shared.EnablePlots(); ProcessControl pc = new ProcessControl(); var dataSet = pc.CascadeControl(); - - Shared.EnablePlots(); + Shared.DisablePlots(); + // Shared.EnablePlots(); } [Test] public void FeedForwardControl_Part1() { - Shared.DisablePlots(); - + // Shared.EnablePlots(); + // Shared.EnablePlots(); ProcessControl pc = new ProcessControl(); var dataSet = pc.FeedForward_Part1(); + // Shared.DisablePlots(); - Shared.EnablePlots(); + Assert.IsTrue(dataSet.GetValue("Process1-Output_Y", 599) -60< 0.01); + + // Shared.EnablePlots(); } diff --git a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/BasicPIDandSisoTests.cs similarity index 88% rename from TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs rename to TimeSeriesAnalysis.Tests/Test/PlantSimulations/BasicPIDandSisoTests.cs index f30f55b3..83b4fae0 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/BasicPIDandSisoTests.cs @@ -14,7 +14,7 @@ namespace TimeSeriesAnalysis.Test.PlantSimulations { [TestFixture] - class SISOTests + class BasicPIDandSISOTests { int timeBase_s = 1; int N = 500; @@ -80,13 +80,14 @@ public void SimulateSingle_InitsRunsAndConverges() inputData.CreateTimestamps(timeBase_s); var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY = simData.GetValues(processModel1.GetID(), SignalType.Output_Y); Assert.IsTrue(Math.Abs(simY[0] - 55) < 0.01); Assert.IsTrue(Math.Abs(simY.Last() - 60) < 0.01); // now test that "simulateSingle" produces the same result! - var isOk2 = plantSim.SimulateSingle(inputData, processModel1.ID, out TimeSeriesDataSet simData2); + var isOk2 = PlantSimulatorHelper.SimulateSingle(inputData, processModel1, out TimeSeriesDataSet simData2); + // var isOk2 = plantSim.SimulateSingle(inputData, processModel1.ID, out TimeSeriesDataSet simData2); if (false) { @@ -193,7 +194,7 @@ public void SimulateSingle_SecondOrderSystem() timeBase_s, TestContext.CurrentContext.Test.Name+ v1); Shared.DisablePlots(); } - CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); } @@ -219,12 +220,13 @@ public void SimulateSingle_NullGains_RunsWithZeroOutput() inputData.CreateTimestamps(timeBase_s); var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY = simData.GetValues(processModel1.GetID(), SignalType.Output_Y); Assert.IsTrue(Math.Abs(simY[0]) == 0); // now test that "simulateSingle" produces the same result! - var isOk2 = plantSim.SimulateSingle(inputData, processModel1.ID, out TimeSeriesDataSet simData2); + var isOk2 = PlantSimulatorHelper.SimulateSingle(inputData, processModel1, out TimeSeriesDataSet simData2); + //var isOk2 = plantSim.SimulateSingle(inputData, processModel1.ID, out TimeSeriesDataSet simData2); //plots @@ -253,110 +255,12 @@ public void SimulateSingle_NullGains_RunsWithZeroOutput() } - [TestCase] - public void Serial2_RunsAndConverges() - { - var plantSim = new PlantSimulator( - new List { processModel1, processModel2 },"Serial2"); - - plantSim.ConnectModels(processModel1, processModel2); - var inputData = new TimeSeriesDataSet(); - inputData.Add(plantSim.AddExternalSignal(processModel1,SignalType.External_U), TimeSeriesCreator.Step(N / 4, N, 50, 55)); - inputData.CreateTimestamps(timeBase_s); - var isOk = plantSim.Simulate(inputData,out TimeSeriesDataSet simData); - plantSim.Serialize(); - Assert.IsTrue(isOk); - CommonAsserts(inputData, simData, plantSim); - double[] simY = simData.GetValues(processModel2.GetID(), SignalType.Output_Y); - Assert.IsTrue(Math.Abs(simY[0] - (55 * 1.1 + 5)) < 0.01); - Assert.IsTrue(Math.Abs(simY.Last() - (60 * 1.1 + 5)) < 0.01); - /* Plot.FromList(new List { - simData.GetValues(processModel1.GetID(),SignalType.Output_Y_sim), - simData.GetValues(processModel2.GetID(),SignalType.Output_Y_sim), - simData.GetValues(processModel1.GetID(),SignalType.External_U)}, - new List { "y1=y_sim1", "y1=y_sim2", "y3=u" }, - timeBase_s, "UnitTest_SerialProcess");*/ - } - [TestCase] - public void Serial3_RunsAndConverges() - { - var plantSim = new PlantSimulator( - new List { processModel1, processModel2, processModel3 }); - - plantSim.ConnectModels(processModel1, processModel2); - plantSim.ConnectModels(processModel2, processModel3); - - var inputData = new TimeSeriesDataSet(); - inputData.Add(plantSim.AddExternalSignal(processModel1, SignalType.External_U), TimeSeriesCreator.Step(N / 4, N, 50, 55)); - inputData.CreateTimestamps(timeBase_s); - var isOk = plantSim.Simulate(inputData,out TimeSeriesDataSet simData); - - //var serialIsOk = plantSim.Serialize("SISO_Serial3",@"c:\appl"); - //Assert.IsTrue(serialIsOk); - - Assert.IsTrue(isOk); - CommonAsserts(inputData, simData, plantSim); - - double[] simY = simData.GetValues(processModel3.GetID(), SignalType.Output_Y); - Assert.IsTrue(Math.Abs(simY[0] - ((55 * 1.1 + 5)*1.1+5)) < 0.01); - Assert.IsTrue(Math.Abs(simY.Last() - ((60 * 1.1 + 5)*1.1+5)) < 0.01); - - if (false) - { - Plot.FromList(new List { - simData.GetValues(processModel1.GetID(),SignalType.Output_Y), - simData.GetValues(processModel2.GetID(),SignalType.Output_Y), - simData.GetValues(processModel3.GetID(),SignalType.Output_Y), - inputData.GetValues(processModel1.GetID(),SignalType.External_U)}, - new List { "y1=y_sim1", "y1=y_sim2", "y1=y_sim3", "y3=u" }, - timeBase_s, TestContext.CurrentContext.Test.Name); - } - } - - - - - - - - - - public static void CommonAsserts(TimeSeriesDataSet inputData,TimeSeriesDataSet simData, PlantSimulator plant) - { - Assert.IsNotNull(simData,"simData should not be null"); - - var signalNames = simData.GetSignalNames(); - - foreach (string signalName in signalNames) - { - var signal = simData.GetValues(signalName); - // test that all systems start in steady-state - double firstTwoValuesDiff = Math.Abs(signal.ElementAt(0) - signal.ElementAt(1)); - double lastTwoValuesDiff = Math.Abs(signal.ElementAt(signal.Length - 2) - signal.ElementAt(signal.Length - 1)); - - Assert.AreEqual(signal.Count(), simData.GetLength(),"all signals should be same length as N"); - Assert.IsTrue(firstTwoValuesDiff < 0.01, "system should start up in steady-state"); - Assert.IsTrue(lastTwoValuesDiff < 0.01, "system should end up in steady-state"); - } - - - - Assert.AreEqual(simData.GetLength(), simData.GetTimeStamps().Count(), "number of timestamps should match number of data points in sim"); - Assert.AreEqual(simData.GetTimeStamps().Last(), inputData.GetTimeStamps().Last(),"datasets should end at same timestamp"); - - /* foreach (var modelKeyValuePair in plant.GetModels()) - { - Assert.IsNotNull(simData.GetValues(modelKeyValuePair.Value.GetID(), SignalType.Output_Y),"model output was not simulated"); - }*/ - - } - private void BasicPIDCommonTests(TimeSeriesDataSet simData, PidModel pidModel) { @@ -469,6 +373,75 @@ public void BasicPID_SetpointStep_RunsAndConverges(bool delayPidOutputOneSample) BasicPIDCommonTests(simData, pidModel); } + + [TestCase] + public void Serial2_SISO_RunsAndConverges() + { + var plantSim = new PlantSimulator( + new List { processModel1, processModel2 }, "Serial2"); + + plantSim.ConnectModels(processModel1, processModel2); + var inputData = new TimeSeriesDataSet(); + inputData.Add(plantSim.AddExternalSignal(processModel1, SignalType.External_U), TimeSeriesCreator.Step(N / 4, N, 50, 55)); + inputData.CreateTimestamps(timeBase_s); + var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); + + plantSim.Serialize(); + + Assert.IsTrue(isOk); + PsTest.CommonAsserts(inputData, simData, plantSim); + + double[] simY = simData.GetValues(processModel2.GetID(), SignalType.Output_Y); + Assert.IsTrue(Math.Abs(simY[0] - (55 * 1.1 + 5)) < 0.01); + Assert.IsTrue(Math.Abs(simY.Last() - (60 * 1.1 + 5)) < 0.01); + + /* Plot.FromList(new List { + simData.GetValues(processModel1.GetID(),SignalType.Output_Y_sim), + simData.GetValues(processModel2.GetID(),SignalType.Output_Y_sim), + simData.GetValues(processModel1.GetID(),SignalType.External_U)}, + new List { "y1=y_sim1", "y1=y_sim2", "y3=u" }, + timeBase_s, "UnitTest_SerialProcess");*/ + } + + + [TestCase] + public void Serial3_SISO_RunsAndConverges() + { + var plantSim = new PlantSimulator( + new List { processModel1, processModel2, processModel3 }); + + plantSim.ConnectModels(processModel1, processModel2); + plantSim.ConnectModels(processModel2, processModel3); + + var inputData = new TimeSeriesDataSet(); + inputData.Add(plantSim.AddExternalSignal(processModel1, SignalType.External_U), TimeSeriesCreator.Step(N / 4, N, 50, 55)); + inputData.CreateTimestamps(timeBase_s); + var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); + + //var serialIsOk = plantSim.Serialize("SISO_Serial3",@"c:\appl"); + //Assert.IsTrue(serialIsOk); + + Assert.IsTrue(isOk); + PsTest.CommonAsserts(inputData, simData, plantSim); + + double[] simY = simData.GetValues(processModel3.GetID(), SignalType.Output_Y); + Assert.IsTrue(Math.Abs(simY[0] - ((55 * 1.1 + 5) * 1.1 + 5)) < 0.01); + Assert.IsTrue(Math.Abs(simY.Last() - ((60 * 1.1 + 5) * 1.1 + 5)) < 0.01); + + if (false) + { + Plot.FromList(new List { + simData.GetValues(processModel1.GetID(),SignalType.Output_Y), + simData.GetValues(processModel2.GetID(),SignalType.Output_Y), + simData.GetValues(processModel3.GetID(),SignalType.Output_Y), + inputData.GetValues(processModel1.GetID(),SignalType.External_U)}, + new List { "y1=y_sim1", "y1=y_sim2", "y1=y_sim3", "y3=u" }, + timeBase_s, TestContext.CurrentContext.Test.Name); + } + } + + + // // Incredibly important that this unit tests passes, as SimulateSingle is used to estimate the disturbance as part of initalization of // Simulate(), so these two methods need to simulate in a consisten way for disturbance estimation to work, which again is vital for @@ -476,7 +449,7 @@ public void BasicPID_SetpointStep_RunsAndConverges(bool delayPidOutputOneSample) // [TestCase] - public void BasicPID_CompareSimulateAndSimulateSingle_MustGiveSameResultForDisturbanceEstToWork() + public void BasicPIDSetpointStep_CompareSimulateAndSimulateSingle_MustGiveSameResultForDisturbanceEstToWork() { //var pidCopy = pidModel1.Clone(); double newSetpoint = 51; @@ -496,16 +469,9 @@ public void BasicPID_CompareSimulateAndSimulateSingle_MustGiveSameResultForDistu newSet.AddSet(simData); newSet.SetTimeStamps(inputData.GetTimeStamps().ToList()); - // slight deviation between SimulateSingle and Simulate regardless of which version is used: - //v1 - // var isOK2 = plantSim.SimulateSingle(newSet, pidModel1.ID,out TimeSeriesDataSet simData2); - //v2 - // pidCopy.SetInputIDs(pidModel1.GetModelInputIDs()); - // var isOk2 = PlantSimulator.SimulateSingle(newSet, pidCopy, out var simData2); - //v3 - var isOk2 = PlantSimulator.SimulateSingle(newSet, pidModel1, out var simData2); + var isOk2 = PlantSimulatorHelper.SimulateSingle(newSet, pidModel1, out var simData2); - if (true) + if (false) { Shared.EnablePlots(); Plot.FromList(new List { @@ -572,5 +538,43 @@ public void BasicPID_SetpointStep_WithNoiseAndFiltering_FilteringWorks() } } + + [TestCase(0)] + [TestCase(1)] + [TestCase(10)] + public void TimeDelay(int timeDelay_s) + { + var timeBase_s = 1; + var parameters = new UnitParameters + { + LinearGains = new double[] { 1 }, + TimeConstant_s = 0, + TimeDelay_s = timeDelay_s, + Bias = 0 + }; + var model = new UnitModel(parameters); + double[] u1 = Vec.Concat(Vec.Fill(0, 31), + Vec.Fill(1, 30)); + UnitDataSet dataSet = new UnitDataSet(); + dataSet.U = Array2D.CreateFromList(new List { u1 }); + dataSet.CreateTimeStamps(timeBase_s); + + (bool isOk, double[] y_sim) = PlantSimulatorHelper.SimulateSingle(dataSet, model); + // plot + bool doPlot = false; + if (doPlot) + { + Shared.EnablePlots(); + Plot.FromList(new List { y_sim, u1 }, new List { "y1=ymeas ", "y3=u1" }, timeBase_s); + Shared.DisablePlots(); + } + //assert + // Assert.IsNotNull(retSim); + Assert.IsTrue(isOk); + Assert.IsTrue(y_sim[30 + timeDelay_s] == 0, "step should not arrive at y_sim too early"); + Assert.IsTrue(y_sim[31 + timeDelay_s] == 1, "steps should be delayed exactly timeDelay_s later "); + } + + } } diff --git a/TimeSeriesAnalysis.Tests/Tests/GainSchedModelTests.cs b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/GainSchedSimulateTests.cs similarity index 98% rename from TimeSeriesAnalysis.Tests/Tests/GainSchedModelTests.cs rename to TimeSeriesAnalysis.Tests/Test/PlantSimulations/GainSchedSimulateTests.cs index 52437ec4..81da0150 100644 --- a/TimeSeriesAnalysis.Tests/Tests/GainSchedModelTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/GainSchedSimulateTests.cs @@ -8,10 +8,10 @@ namespace TimeSeriesAnalysis.Test.PlantSimulations { /// - /// Test of process simulations where each of or some of the models have multiple inputs + /// Test of PlantSimulator simulations of the GainSched model /// [TestFixture] - class GainsSchedModelTests + class GainsSchedSimulateTests { int timeBase_s = 1; int N = 480; @@ -43,7 +43,7 @@ class GainsSchedModelTests TimeDelay_s = 1, GainSchedParameterIndex = 0 }; - GainSchedParameters gainSchedP4_singleThreshold_singleInput_bias_and_timedelay= new GainSchedParameters(0,-15) + GainSchedParameters gainSchedP4_singleThreshold_singleInput_bias_and_timedelay = new GainSchedParameters(0,-15) { TimeConstant_s = new double[] { 10, 20 }, TimeConstantThresholds = new double[] { 2 }, @@ -170,7 +170,7 @@ public void TenThresholds_DifferentGainsAboveEachThreshold(bool useHighOperating // Act var isSimulatable = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY1 = simData.GetValues(refModel.GetID(), SignalType.Output_Y); if (false) @@ -316,7 +316,7 @@ public void ContinousGradualRamp_BumplessModelOutput(int ver, string upOrDown) timeBase_s, "ContinousGradualRamp_BumplessModelOutput" + ver + upOrDown); Shared.DisablePlots(); } - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); // assert that step increase is gradual, even at the boundaries between different double maxChg = ((11.0 +1)* 10.0) / N; @@ -361,8 +361,8 @@ public void SingleThreshold_DifferentGainsAboveAndBelowThreshold(int ver, double // Act var isSimulatable = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); - - SISOTests.CommonAsserts(inputData, simData, plantSim); + + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY1 = simData.GetValues(gainSched.GetID(), SignalType.Output_Y); bool doPlot = false;// should be false unless debugging. @@ -433,7 +433,7 @@ public void NonzeroOperatingPoint_SimulationStartsInOpPointOk(double uOperatingP Shared.DisablePlots(); } - SISOTests.CommonAsserts(inputData, refSimData, truePlantSim); + PsTest.CommonAsserts(inputData, refSimData, truePlantSim); Assert.That(unitData.Y_meas.First() == yOperatingPoint, "since time series starts in uOperatingPoint, simulation should start in yOperatingPoint"); } diff --git a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorModelTests.cs b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/LargerSystemSimulations.cs similarity index 94% rename from TimeSeriesAnalysis.Tests/Tests/PlantSimulatorModelTests.cs rename to TimeSeriesAnalysis.Tests/Test/PlantSimulations/LargerSystemSimulations.cs index 928fbc9e..5cde022a 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorModelTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/LargerSystemSimulations.cs @@ -11,7 +11,7 @@ namespace TimeSeriesAnalysis.Test.PlantSimulations /// Test of process simulations where each of or some of the models have multiple inputs /// [TestFixture] - class PlantSimulatorModelTests + class LargerSystemSimulations { int timeBase_s = 1; int N = 480; @@ -111,7 +111,7 @@ public void MinSelectWithPID_RunsAndConverges() var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY = simData.GetValues(minSelect1.GetID(), SignalType.SelectorOut); SerializeHelper.Serialize("MinSelectWithPID", plantSim, inputData, simData); @@ -141,7 +141,7 @@ public void MaxSelect_RunsAndConverges() var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY = simData.GetValues(maxSelect1.GetID(), SignalType.SelectorOut); Assert.IsTrue(Math.Abs(simY[0] - (6.7)) < 0.01); @@ -164,7 +164,7 @@ public void Single_RunsAndConverges() inputData.CreateTimestamps(timeBase_s); var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY = simData.GetValues(processModel1.GetID(), SignalType.Output_Y); Assert.IsTrue(Math.Abs(simY[0] - (1 * 50 + 0.5 * 50 + 5)) < 0.01); @@ -221,7 +221,7 @@ public void PIDAndSingle_RunsAndConverges(bool doReverseInputConnections) } */ [TestCase] - public void Serial2_RunsAndConverges() + public void Serial2_MISO_RunsAndConverges() { var plantSim = new PlantSimulator( new List { processModel1, processModel2 }); @@ -236,7 +236,7 @@ public void Serial2_RunsAndConverges() var isOk = plantSim.Simulate(inputData,out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData,simData, plantSim); + PsTest.CommonAsserts(inputData,simData, plantSim); //Shared.EnablePlots(); //Plot.FromList(new List { @@ -311,18 +311,18 @@ public void PIDandSerial2_RunsAndConverges(int ver) } Assert.IsTrue(isOk,"simulation returned false, it failed"); - /* - Plot.FromList(new List { - simData.GetValues(processModel1.GetID(),SignalType.Output_Y_sim), - simData.GetValues(processModel2.GetID(),SignalType.Output_Y_sim), - simData.GetValues(pidModel1.GetID(),SignalType.PID_U), - simData.GetValues(processModel1.GetID(),SignalType.External_U,externalUIndex), - simData.GetValues(processModel2.GetID(),SignalType.External_U,(int)INDEX.SECOND) - }, - new List { "y1=y_sim1","y1=y_sim2", "y3=u1(pid)", "y3=u2", "y3=u3" }, - timeBase_s, "UnitTest_PIDandSerial2");*/ - - SISOTests.CommonAsserts(inputData, simData, plantSim); + /* + Plot.FromList(new List { + simData.GetValues(processModel1.GetID(),SignalType.Output_Y_sim), + simData.GetValues(processModel2.GetID(),SignalType.Output_Y_sim), + simData.GetValues(pidModel1.GetID(),SignalType.PID_U), + simData.GetValues(processModel1.GetID(),SignalType.External_U,externalUIndex), + simData.GetValues(processModel2.GetID(),SignalType.External_U,(int)INDEX.SECOND) + }, + new List { "y1=y_sim1","y1=y_sim2", "y3=u1(pid)", "y3=u2", "y3=u3" }, + timeBase_s, "UnitTest_PIDandSerial2");*/ + + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY = simData.GetValues(processModel2.GetID(), SignalType.Output_Y); Assert.IsTrue(Math.Abs(simY[0] - 150) < 0.01, "unexpected starting value"); Assert.IsTrue(Math.Abs(simY.Last() - 150) < 0.1, "unexpected ending value"); @@ -357,7 +357,7 @@ public void ComputationalLoop_TwoModels_RunsAndConverges() var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); /* Shared.EnablePlots(); Plot.FromList(new List { @@ -394,7 +394,7 @@ public void ComputationalLoop_TwoModelsLoop_One_Upstream_RunsAndConverges() var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); /* Shared.EnablePlots(); Plot.FromList(new List { @@ -434,7 +434,7 @@ public void ComputationalLoop_TwoModelsLoop_OneUpstream_OneDownstream_RunsAndCon var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); /* Shared.EnablePlots(); @@ -475,7 +475,7 @@ public void ComputationalLoop_ThreeModelsLoop_RunsAndConverges() var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); /* Shared.EnablePlots(); Plot.FromList(new List { @@ -490,13 +490,18 @@ public void ComputationalLoop_ThreeModelsLoop_RunsAndConverges() Shared.DisablePlots();*/ } + + + + + [TestCase(0,Description ="this is the easiest, as it requires no res-")] [TestCase(1)] [TestCase(2)] [TestCase(2)] - public void Serial3_RunsAndConverges(int ver) + public void Serial3_MISO_RunsAndConverges(int ver) { List modelList = new List(); if (ver == 0) @@ -533,7 +538,7 @@ public void Serial3_RunsAndConverges(int ver) var isOk = plantSim.Simulate(inputData,out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); double[] simY = simData.GetValues(processModel3.GetID(), SignalType.Output_Y); double expStartVal = ((1 * 50 + 0.5 * 50 + 5) * 1.1 + 50 * 0.6 + 5)*0.8 + 0.7*30 + 5; @@ -583,7 +588,7 @@ public void Divide_RunsAndConverges() var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); SerializeHelper.Serialize("Divide", plantSim, inputData, simData); @@ -611,7 +616,7 @@ public void TwoProcessesToOne_RunsAndConverges() var isOk = plantSim.Simulate(inputData, out TimeSeriesDataSet simData); Assert.IsTrue(isOk); - SISOTests.CommonAsserts(inputData, simData, plantSim); + PsTest.CommonAsserts(inputData, simData, plantSim); SerializeHelper.Serialize("TwoProcessesToOne",plantSim,inputData,simData); diff --git a/TimeSeriesAnalysis.Tests/Test/PlantSimulations/PsTest.cs b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/PsTest.cs new file mode 100644 index 00000000..3c07ca3d --- /dev/null +++ b/TimeSeriesAnalysis.Tests/Test/PlantSimulations/PsTest.cs @@ -0,0 +1,48 @@ +using NUnit.Framework; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using TimeSeriesAnalysis.Dynamic; + +namespace TimeSeriesAnalysis.Test.PlantSimulations +{ + /// + /// Common helper class for PlantSimulator tests + /// + internal class PsTest + { + + + public static void CommonAsserts(TimeSeriesDataSet inputData, TimeSeriesDataSet simData, PlantSimulator plant) + { + Assert.IsNotNull(simData, "simData should not be null"); + + var signalNames = simData.GetSignalNames(); + + foreach (string signalName in signalNames) + { + var signal = simData.GetValues(signalName); + // test that all systems start in steady-state + double firstTwoValuesDiff = Math.Abs(signal.ElementAt(0) - signal.ElementAt(1)); + double lastTwoValuesDiff = Math.Abs(signal.ElementAt(signal.Length - 2) - signal.ElementAt(signal.Length - 1)); + + Assert.AreEqual(signal.Count(), simData.GetLength(), "all signals should be same length as N"); + Assert.IsTrue(firstTwoValuesDiff < 0.01, "system should start up in steady-state"); + Assert.IsTrue(lastTwoValuesDiff < 0.01, "system should end up in steady-state"); + } + + Assert.AreEqual(simData.GetLength(), simData.GetTimeStamps().Count(), "number of timestamps should match number of data points in sim"); + Assert.AreEqual(simData.GetTimeStamps().Last(), inputData.GetTimeStamps().Last(), "datasets should end at same timestamp"); + + /* foreach (var modelKeyValuePair in plant.GetModels()) + { + Assert.IsNotNull(simData.GetValues(modelKeyValuePair.Value.GetID(), SignalType.Output_Y),"model output was not simulated"); + }*/ + + } + + + } +} diff --git a/TimeSeriesAnalysis.Tests/Tests/PlotUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/PlotUnitTests.cs similarity index 100% rename from TimeSeriesAnalysis.Tests/Tests/PlotUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/PlotUnitTests.cs diff --git a/TimeSeriesAnalysis.Tests/Tests/CsvTest.cs b/TimeSeriesAnalysis.Tests/Test/Serialization/CsvTest.cs similarity index 98% rename from TimeSeriesAnalysis.Tests/Tests/CsvTest.cs rename to TimeSeriesAnalysis.Tests/Test/Serialization/CsvTest.cs index f809ee79..6a3ef31f 100644 --- a/TimeSeriesAnalysis.Tests/Tests/CsvTest.cs +++ b/TimeSeriesAnalysis.Tests/Test/Serialization/CsvTest.cs @@ -8,7 +8,7 @@ using NUnit.Framework; using NUnit.Framework.Interfaces; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.Serialization { [TestFixture] class FileReading diff --git a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSerialization.cs b/TimeSeriesAnalysis.Tests/Test/Serialization/PlantSimulatorSerialization.cs similarity index 99% rename from TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSerialization.cs rename to TimeSeriesAnalysis.Tests/Test/Serialization/PlantSimulatorSerialization.cs index b0e9062d..90e70b2e 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSerialization.cs +++ b/TimeSeriesAnalysis.Tests/Test/Serialization/PlantSimulatorSerialization.cs @@ -8,7 +8,7 @@ using TimeSeriesAnalysis.Dynamic; using TimeSeriesAnalysis.Utility; -namespace TimeSeriesAnalysis.Tests.Serialization +namespace TimeSeriesAnalysis.Test.Serialization { internal class PlantSimulatorSerialization { diff --git a/TimeSeriesAnalysis.Tests/Tests/CorrelationCalculator.cs b/TimeSeriesAnalysis.Tests/Test/SysID/CorrelationCalculator.cs similarity index 99% rename from TimeSeriesAnalysis.Tests/Tests/CorrelationCalculator.cs rename to TimeSeriesAnalysis.Tests/Test/SysID/CorrelationCalculator.cs index dae10528..4c9a64ad 100644 --- a/TimeSeriesAnalysis.Tests/Tests/CorrelationCalculator.cs +++ b/TimeSeriesAnalysis.Tests/Test/SysID/CorrelationCalculator.cs @@ -8,7 +8,7 @@ using TimeSeriesAnalysis.Utility; using NUnit.Framework; -namespace TimeSeriesAnalysis.Test +namespace TimeSeriesAnalysis.Test.SysID { [TestFixture] class CorrelationCalculatorTests diff --git a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs b/TimeSeriesAnalysis.Tests/Test/SysID/GainSchedIdentifyTests.cs similarity index 96% rename from TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs rename to TimeSeriesAnalysis.Tests/Test/SysID/GainSchedIdentifyTests.cs index 73ff02c3..2349121a 100644 --- a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/SysID/GainSchedIdentifyTests.cs @@ -158,7 +158,7 @@ public void FiveGainsStatic_StepChange_ForGivenThresholds_CorrectGains(int ver, var dataSet = new UnitDataSet(); dataSet.SetU(input); dataSet.CreateTimeStamps(timeBase_s); - PlantSimulator.SimulateSingleToYmeas(dataSet, refModel,noiseAmplitude); + PlantSimulatorHelper.SimulateSingleToYmeas(dataSet, refModel,noiseAmplitude); var idModel = GainSchedIdentifier.IdentifyForGivenThresholds(dataSet, gsFittingSpecs); @@ -243,7 +243,7 @@ public void IgnoreIndicesInMiddleOfDataset_ResultShouldStillBeGood( double fitSc unitData1.CreateTimeStamps(timeBase_s); GainSchedModel trueModel = new GainSchedModel(trueGSparams, "true gain sched model"); - PlantSimulator.SimulateSingleToYmeas(unitData1, trueModel, noiseAmp, 454); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData1, trueModel, noiseAmp, 454); } // Dataset 2 var unitData2 = new UnitDataSet("dataset2"); @@ -260,7 +260,7 @@ public void IgnoreIndicesInMiddleOfDataset_ResultShouldStillBeGood( double fitSc TimeDelay_s = timeBase_s * 0, }; GainSchedModel trueModel = new GainSchedModel(otherGsParams, "true gain sched model"); - PlantSimulator.SimulateSingleToYmeas(unitData2, trueModel, noiseAmp, 454); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData2, trueModel, noiseAmp, 454); } // dataset 3 var unitData3 = new UnitDataSet("dataset3"); @@ -269,7 +269,7 @@ public void IgnoreIndicesInMiddleOfDataset_ResultShouldStillBeGood( double fitSc unitData3.SetU(u3); unitData3.CreateTimeStamps(timeBase_s); GainSchedModel trueModel = new GainSchedModel(trueGSparams, "true gain sched model"); - PlantSimulator.SimulateSingleToYmeas(unitData3, trueModel, noiseAmp, 454); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData3, trueModel, noiseAmp, 454); } var joinedDataSet = new UnitDataSet(unitData1); @@ -335,7 +335,7 @@ public void TimeDelaySingleTc_StepChange_Identify_TcAndTdEstOk(int timeDelaySamp }; GainSchedModel trueModel = new GainSchedModel(trueGSparams, "true gain sched model"); - PlantSimulator.SimulateSingleToYmeas(unitData, trueModel, noiseAmp, 454); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData, trueModel, noiseAmp, 454); // Act var idModel = GainSchedIdentifier.Identify(unitData); @@ -395,7 +395,7 @@ public void NonzeroOperatingPointU_Both_EstimatesStillOk(double uOperatingPoint, }; GainSchedModel trueModel = new GainSchedModel(trueGSparams, "True Model"); - PlantSimulator.SimulateSingleToYmeas(unitData, trueModel, noiseAmp, (int)Math.Ceiling(2 * gainSchedThreshold + 45)); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData, trueModel, noiseAmp, (int)Math.Ceiling(2 * gainSchedThreshold + 45)); // Act var idModel = new GainSchedModel(); @@ -458,7 +458,7 @@ public void TwoGainsConstTc_RampChange_Both_AllParamsEstOk(string solver, double var unitData = new UnitDataSet(); unitData.SetU(input); unitData.CreateTimeStamps(timeBase_s); - bool isOk = PlantSimulator.SimulateSingleToYmeas(unitData,trueModel, noise_abs); + bool isOk = PlantSimulatorHelper.SimulateSingleToYmeas(unitData,trueModel, noise_abs); GainSchedModel idModel = new GainSchedModel(); if (solver == "Identify") @@ -542,7 +542,7 @@ public void TwoGainsAndTwoTc_StepChange_Identify_TwoTcEstEstOk(double gain_sched }; GainSchedModel trueModel = new GainSchedModel(trueGSparams, "Correct gain sched model"); - PlantSimulator.SimulateSingleToYmeas(unitData,trueModel, noiseAmp, (int)Math.Ceiling(2 * gain_sched_threshold + 45)); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData,trueModel, noiseAmp, (int)Math.Ceiling(2 * gain_sched_threshold + 45)); // Act var idModel = GainSchedIdentifier.Identify(unitData); @@ -606,13 +606,13 @@ public void ChangeOperatingPoint_YsimUnchanged(double origOpU, double origOpY) // make the bias nonzero to test that the operating point estimation works. GainSchedModel trueModel = new GainSchedModel(trueParams, "true"); - PlantSimulator.SimulateSingleToYsim(unitData, trueModel); + PlantSimulatorHelper.SimulateSingleToYsim(unitData, trueModel); var alteredParams = new GainSchedParameters(trueModel.GetModelParameters()); var alteredIdModel = new GainSchedModel(alteredParams,"altered"); alteredIdModel.GetModelParameters().MoveOperatingPointUWithoutChangingModel(3); - (bool isOk3, double[] simY2) = PlantSimulator.SimulateSingle(unitData, alteredIdModel); + (bool isOk3, double[] simY2) = PlantSimulatorHelper.SimulateSingle(unitData, alteredIdModel); // plot bool doPlot = false; @@ -668,7 +668,7 @@ public void FlatData_IdTerminatesWithoutCrashing(string solverId) // trueParams.OperatingPoint_Y = 1.34; GainSchedModel trueModel = new GainSchedModel(trueParams, "Correct gain sched model"); - PlantSimulator.SimulateSingleToYmeas(unitData, trueModel, noiseAmplitude, 123); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData, trueModel, noiseAmplitude, 123); // Act var idModel = new GainSchedModel(); @@ -738,9 +738,9 @@ public void TwoGainsConstTc_StepChange_Both_TCAndThresholdFoundOk(double TimeCon // make the bias nonzero to test that the operating point estimation works. // trueParams.OperatingPoint_Y = 1.34; - GainSchedModel trueModel = new GainSchedModel(trueParams, "Correct gain sched model"); + var trueModel = new GainSchedModel(trueParams, "Correct gain sched model"); - PlantSimulator.SimulateSingleToYmeas(unitData, trueModel, noiseAmplitude,123); + PlantSimulatorHelper.SimulateSingleToYmeas(unitData, trueModel, noiseAmplitude,123); // Act var idModel = new GainSchedModel(); @@ -756,11 +756,11 @@ public void TwoGainsConstTc_StepChange_Both_TCAndThresholdFoundOk(double TimeCon gsFittingSpecs.uTimeConstantThresholds = trueParams.TimeConstantThresholds; idModel = GainSchedIdentifier.IdentifyForGivenThresholds(unitData, gsFittingSpecs); } - (bool isOk2, double[] simY2) = PlantSimulator.SimulateSingle(unitData, idModel); + (bool isOk2, double[] simY2) = PlantSimulatorHelper.SimulateSingle(unitData, idModel); var alteredIdModel = (GainSchedModel)idModel.Clone("altered"); alteredIdModel.GetModelParameters().MoveOperatingPointUWithoutChangingModel(3); - (bool isOk3, double[] simY3) = PlantSimulator.SimulateSingle(unitData, alteredIdModel); + (bool isOk3, double[] simY3) = PlantSimulatorHelper.SimulateSingle(unitData, alteredIdModel); // plot bool doPlot = false; diff --git a/TimeSeriesAnalysis.Tests/Tests/PidIdentUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/SysID/PidIdentUnitTests.cs similarity index 99% rename from TimeSeriesAnalysis.Tests/Tests/PidIdentUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/SysID/PidIdentUnitTests.cs index 35831f99..74c4fa32 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PidIdentUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/SysID/PidIdentUnitTests.cs @@ -5,7 +5,7 @@ using TimeSeriesAnalysis.Dynamic; using TimeSeriesAnalysis.Utility; -namespace TimeSeriesAnalysis.Test.PidID +namespace TimeSeriesAnalysis.Test.SysID { [TestFixture] class PidIdentUnitTests diff --git a/TimeSeriesAnalysis.Tests/Tests/SysIDUnitTests.cs b/TimeSeriesAnalysis.Tests/Test/SysID/UnitIdentificiationTests.cs similarity index 96% rename from TimeSeriesAnalysis.Tests/Tests/SysIDUnitTests.cs rename to TimeSeriesAnalysis.Tests/Test/SysID/UnitIdentificiationTests.cs index 4efb72a2..30f28a68 100644 --- a/TimeSeriesAnalysis.Tests/Tests/SysIDUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Test/SysID/UnitIdentificiationTests.cs @@ -13,55 +13,13 @@ using System.Xml; namespace TimeSeriesAnalysis.Test.SysID -{ - - class UnitSimulation - { - [TestCase(0)] - [TestCase(1)] - [TestCase(10)] - public void UnitSimulate_TimeDelay(int timeDelay_s) - { - var timeBase_s = 1; - var parameters = new UnitParameters - { - LinearGains = new double []{ 1}, - TimeConstant_s = 0, - TimeDelay_s = timeDelay_s, - Bias =0 - }; - var model = new UnitModel(parameters); - double[] u1 = Vec.Concat(Vec.Fill(0, 31), - Vec.Fill(1, 30)); - UnitDataSet dataSet = new UnitDataSet(); - dataSet.U = Array2D.CreateFromList(new List { u1 }); - dataSet.CreateTimeStamps(timeBase_s); - - (bool isOk,double[] y_sim) = PlantSimulator.SimulateSingle(dataSet, model); - // plot - bool doPlot = false; - if (doPlot) - { - Shared.EnablePlots(); - Plot.FromList(new List { y_sim, u1 }, new List { "y1=ymeas ", "y3=u1" }, timeBase_s); - Shared.DisablePlots(); - } - //assert - // Assert.IsNotNull(retSim); - Assert.IsTrue(isOk); - Assert.IsTrue(y_sim[30+ timeDelay_s] == 0,"step should not arrive at y_sim too early"); - Assert.IsTrue(y_sim[31+ timeDelay_s] == 1, "steps should be delayed exactly timeDelay_s later "); - } - - } - - +{ /// /// DefaultModel unit tests /// In the naming convention I1 refers to one input, I2 two inputs etc. /// - class UnitIdentification + class UnitIdentificationTests { static Plot4Test plot = new Plot4Test(false); double timeBase_s = 1; @@ -71,18 +29,16 @@ public UnitDataSet CreateDataSet(UnitParameters designParameters, double[,] U, d double noiseAmplitude = 0, bool addInBadDataToYmeasAndU = false, double badValueId = Double.NaN, bool doNonWhiteNoise=false) { - UnitModel model = new UnitModel(designParameters); + var model = new UnitModel(designParameters); this.timeBase_s = timeBase_s; - UnitDataSet dataSet = new UnitDataSet(); + var dataSet = new UnitDataSet(); dataSet.U = U; dataSet.BadDataID = badValueId; dataSet.CreateTimeStamps(timeBase_s); - // var simulator = new UnitSimulator(model); if (doNonWhiteNoise) { - PlantSimulator.SimulateSingleToYmeas(dataSet, model, 0); - // simulator.SimulateYmeas(ref dataSet, 0); + PlantSimulatorHelper.SimulateSingleToYmeas(dataSet, model, 0); double rand = 0; var randObj = new Random(45466545); for (int i = 0; i < dataSet.Y_meas.Length; i++) @@ -93,8 +49,7 @@ public UnitDataSet CreateDataSet(UnitParameters designParameters, double[,] U, d } else { - PlantSimulator.SimulateSingleToYmeas(dataSet, model, noiseAmplitude); - // simulator.SimulateYmeas(ref dataSet, noiseAmplitude); + PlantSimulatorHelper.SimulateSingleToYmeas(dataSet, model, noiseAmplitude); } if (addInBadDataToYmeasAndU) diff --git a/TimeSeriesAnalysis.Tests/Tests/TimeSeriesDataSetTest.cs b/TimeSeriesAnalysis.Tests/Test/TimeSeriesData/TimeSeriesDataSetTest.cs similarity index 100% rename from TimeSeriesAnalysis.Tests/Tests/TimeSeriesDataSetTest.cs rename to TimeSeriesAnalysis.Tests/Test/TimeSeriesData/TimeSeriesDataSetTest.cs diff --git a/TimeSeriesAnalysis.Tests/Tests/UnitSimulatorTests.cs b/TimeSeriesAnalysis.Tests/Test/UnitSimulatorTests.cs similarity index 100% rename from TimeSeriesAnalysis.Tests/Tests/UnitSimulatorTests.cs rename to TimeSeriesAnalysis.Tests/Test/UnitSimulatorTests.cs diff --git a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSingleTests.cs b/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSingleTests.cs deleted file mode 100644 index 7594549c..00000000 --- a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSingleTests.cs +++ /dev/null @@ -1,218 +0,0 @@ -using NUnit.Framework; -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using System.Threading.Tasks; -using TimeSeriesAnalysis.Dynamic; -using TimeSeriesAnalysis.Utility; - -namespace TimeSeriesAnalysis.Test.DisturbanceID -{ - internal class PlantSimulatorSingleTests - { - - UnitParameters staticModelParameters = new UnitParameters - { - TimeConstant_s = 0, - LinearGains = new double[] { 1.5 }, - TimeDelay_s = 0, - Bias = 5 - }; - - UnitParameters dynamicModelParameters = new UnitParameters - { - TimeConstant_s = 10, - LinearGains = new double[] { 1.5 }, - TimeDelay_s = 5, - Bias = 5 - }; - - PidParameters pidParameters1 = new PidParameters() - { - Kp = 0.2, - Ti_s = 20 - }; - - int timeBase_s = 1; - int N = 300; - DateTime t0 = new DateTime(2010, 1, 1); - - // an extension of the above test to use the more general PlantSimulator.Simulate, rather than the PlantSimulator.SimulateSingle - [TestCase(4),Explicit] - public void PlantSimulator_StepDisturbance_EstimatesOk(double stepAmplitude) - { - var locModelParameters = new UnitParameters - { - TimeConstant_s = 15, - LinearGains = new double[] { 1 }, - TimeDelay_s = 0, - Bias = 5 - }; - - var trueDisturbance = TimeSeriesCreator.Step(100, N, 0, stepAmplitude); - DisturbanceTestUsingPlantSimulator(locModelParameters, trueDisturbance); - } - - [TestCase(4, 1),Explicit] - [TestCase(4, -1)] - public void PlantSimulator_StepDisturbanceANDSetPointStep_EstimatesOk(double disturbanceStepAmplitude, double setpointAmplitude) - { - // Shared.EnablePlots(); - var trueDisturbance = TimeSeriesCreator.Step(100, N, 0, disturbanceStepAmplitude); - var setpoint = TimeSeriesCreator.Step(50, N, 50, 50 + setpointAmplitude); - DisturbanceTestUsingPlantSimulator(dynamicModelParameters, trueDisturbance, true, setpoint); - //Shared.DisablePlots(); - } - - [TestCase(-4), Explicit] - [TestCase(4)] - public void PlantSimulatorSingle_StepDisturbance_EstimatesOk(double stepAmplitude) - { - var trueDisturbance = TimeSeriesCreator.Step(100, N, 0, stepAmplitude); - DisturbanceTestUsingPlantSimulateSingle(new UnitModel(dynamicModelParameters, "PlantSim_d"), trueDisturbance); - } - - - - public void DisturbanceTestUsingPlantSimulator(UnitParameters unitParams, double[] trueDisturbance, - bool doAssertResult = true, double[] externalYset = null) - { - TimeSeriesDataSet referenceSimDataSet; - TimeSeriesDataSet referenceInputDataSet; - // 1 .create synthetic dataset - where a "true" known disturbance is specified - { - var pidModel1 = new PidModel(pidParameters1, "PID1"); - var processModel2 = new UnitModel(unitParams, "Proc1"); - var plantSim = new PlantSimulator( - new List { pidModel1, processModel2 }); - plantSim.ConnectModels(processModel2, pidModel1); - plantSim.ConnectModels(pidModel1, processModel2); - var refYsetSignal = "yset"; - plantSim.AddAndConnectExternalSignal(pidModel1, refYsetSignal, SignalType.Setpoint_Yset); - var refDistSignal = "dist"; - plantSim.AddAndConnectExternalSignal(processModel2, refDistSignal, SignalType.Disturbance_D); - - referenceInputDataSet = new TimeSeriesDataSet(); - if (externalYset == null) - referenceInputDataSet.Add(refYsetSignal, TimeSeriesCreator.Constant(50, N)); - else - referenceInputDataSet.Add(refYsetSignal, externalYset); - referenceInputDataSet.Add(refDistSignal, trueDisturbance); - referenceInputDataSet.CreateTimestamps(timeBase_s); - - var simOk = plantSim.Simulate(referenceInputDataSet, out referenceSimDataSet); - Assert.IsTrue(simOk, "simulating reference case failed!"); - } - // 2.create plant model without disturbance, and try to to find the disturbance signal - { - var pidModel1 = new PidModel(pidParameters1, "PID1"); - - var processModel3 = new UnitModel(unitParams, "Proc1"); - - var distSignal = SignalNamer.EstDisturbance(processModel3); - processModel3.AddSignalToOutput(distSignal); - var plantSim = new PlantSimulator( - new List { pidModel1, processModel3 }); - plantSim.ConnectModels(processModel3, pidModel1); - plantSim.ConnectModels(pidModel1, processModel3); - - // signals can really be named anything, but important for this to work that the names are the same - // in the model objects and in the inputData object - var ysetSignal = SignalNamer.GetSignalName(pidModel1.GetID(), SignalType.Setpoint_Yset); - plantSim.AddAndConnectExternalSignal(pidModel1, ysetSignal, SignalType.Setpoint_Yset); - var ymeasSignal = processModel3.GetOutputID(); - plantSim.AddAndConnectExternalSignal(processModel3, ymeasSignal, SignalType.Output_Y); - var uMeasSignal = SignalNamer.GetSignalName(pidModel1.GetID(), SignalType.PID_U); ; - plantSim.AddAndConnectExternalSignal(pidModel1, uMeasSignal, SignalType.PID_U); - // nb! do not specify the disturbance in this case, instead add the "output_Y" from the abo - - ///////////////// - /// - ///adding u and y to inputdata, should enable the plant simulator to back-calculate the disturbance. - /// - var inputData = new TimeSeriesDataSet(); - inputData.Add(ysetSignal, referenceInputDataSet.GetValues("yset")); - inputData.Add(uMeasSignal, referenceSimDataSet.GetValues("PID1", SignalType.PID_U)); - inputData.Add(ymeasSignal, referenceSimDataSet.GetValues("Proc1", SignalType.Output_Y)); - ///////////////// - Assert.IsTrue(inputData.GetSignalNames().Count() == 3, "configuration errors");//sanity check for configuration errors - inputData.CreateTimestamps(timeBase_s); - - ////////////////////////////////// - var isOK = plantSim.Simulate(inputData, out var simDataSetWithDisturbance); - ////////////////////////////////// - - if (doAssertResult) - { - var pidDataSet = plantSim.GetUnitDataSetForPID(inputData.Combine(simDataSetWithDisturbance), pidModel1); - DisturbanceCalculatorTests.CommonPlotAndAsserts(pidDataSet, simDataSetWithDisturbance.GetValues(distSignal), - trueDisturbance); - } - Assert.IsTrue(isOK, "simulate dataset with the disturbance FAILED!!"); - Assert.IsTrue(plantSim.PlantFitScore == 100, "expect plant fit score 100"); - Assert.IsTrue(simDataSetWithDisturbance.ContainsSignal(distSignal), "simulated dataset does not contain the disturbance signal"); - } - } - - public void DisturbanceTestUsingPlantSimulateSingle(UnitModel processModel, double[] trueDisturbance, - bool doAssertResult = true) - { - TimeSeriesDataSet referenceSimDataSet; - TimeSeriesDataSet referenceInputDataSet; - // 1 .create synthetic dataset - where disturbance is specified - { - var pidModel1 = new PidModel(pidParameters1, "PID1"); - var plantSim = new PlantSimulator( - new List { pidModel1, processModel }); - plantSim.ConnectModels(processModel, pidModel1); - plantSim.ConnectModels(pidModel1, processModel); - referenceInputDataSet = new TimeSeriesDataSet(); - referenceInputDataSet.Add(plantSim.AddExternalSignal(pidModel1, SignalType.Setpoint_Yset), TimeSeriesCreator.Constant(50, N)); - referenceInputDataSet.Add(plantSim.AddExternalSignal(processModel, SignalType.Disturbance_D), trueDisturbance); - referenceInputDataSet.CreateTimestamps(timeBase_s); - var isOK = plantSim.Simulate(referenceInputDataSet, out referenceSimDataSet); - Assert.IsTrue(isOK); - } - // 2.create plant model without disturbance, and try to to find the disturbance signal - { - var pidModel1 = new PidModel(pidParameters1, "PID1"); - var plantSim = new PlantSimulator( - new List { pidModel1, processModel }); - plantSim.ConnectModels(processModel, pidModel1); - plantSim.ConnectModels(pidModel1, processModel); - var inputData = new TimeSeriesDataSet(); - inputData.Add(SignalNamer.GetSignalName(pidModel1.GetID(), SignalType.Setpoint_Yset), - referenceInputDataSet.GetValues(pidModel1.ID, SignalType.Setpoint_Yset)); - ///////////////// - /// - ///adding u and y to inputdata, should enable the plant simulator to back-calculate the disturbance. - /// - // u - inputData.Add(SignalNamer.GetSignalName(pidModel1.GetID(), SignalType.PID_U), - referenceSimDataSet.GetValues(pidModel1.ID, SignalType.PID_U));// use the input u from the other dataset, simulating a "field data" set - //y_meas - should trigger determining the disturbance - inputData.Add(SignalNamer.GetSignalName(processModel.GetID(), SignalType.Output_Y), - referenceSimDataSet.GetValues(processModel.ID, SignalType.Output_Y)); - ///////////////// - inputData.CreateTimestamps(timeBase_s); - var isOK = plantSim.SimulateSingle(inputData, processModel.ID, - out TimeSeriesDataSet simDataSetWithDisturbance); - // var isOK = PlantSimulator.SimulateSingle(inputData, processModel, - // out TimeSeriesDataSet simDataSetWithDisturbance); - - Assert.IsTrue(isOK); - Assert.IsTrue(simDataSetWithDisturbance.ContainsSignal(SignalNamer.EstDisturbance(processModel))); - Assert.IsTrue(simDataSetWithDisturbance.ContainsSignal(processModel.GetID()), - "simData should contain internal output of process as well"); - if (doAssertResult) - { - var pidDataSet = plantSim.GetUnitDataSetForPID(inputData.Combine(simDataSetWithDisturbance), pidModel1); - DisturbanceCalculatorTests.CommonPlotAndAsserts(pidDataSet, simDataSetWithDisturbance.GetValues(SignalNamer.EstDisturbance(processModel)), - trueDisturbance); - } - } - } - - } -} diff --git a/TimeSeriesAnalysis/TimeSeriesDataSet.cs b/TimeSeriesAnalysis/TimeSeriesDataSet.cs index 6a8652b8..be3dabaa 100644 --- a/TimeSeriesAnalysis/TimeSeriesDataSet.cs +++ b/TimeSeriesAnalysis/TimeSeriesDataSet.cs @@ -472,7 +472,7 @@ public string[] GetSignalNames() } /// - /// Get a vector of the timestamps of the data-set + /// Get a vector of the timestamps of the data-set, or null if no timestamps /// /// public DateTime[] GetTimeStamps() diff --git a/docs/articles/processsimulator.md b/docs/articles/processsimulator.md index c077f125..bbed61dc 100644 --- a/docs/articles/processsimulator.md +++ b/docs/articles/processsimulator.md @@ -67,6 +67,8 @@ simulate the disturbance vector as part of the initialization of ``PlantSimulato ### Closed loops : simulation order and disturbance +PID-loops are a special case of a computational loop. + In a closed loop the simulation order will be - *PidModel* reads ``y_meas[k]`` and ``y_setpoint[k]`` and calculates ``u_pid[k]`` @@ -87,4 +89,12 @@ are implemented like this, but often the analysis is done on down-sampled data, **Implicitly the above also defines how to interpret the disturbance d[k].** To be extremely precise with how this is defined is important, as the PlantSimulator is used internally to back-calculte disturbances as is described in the above section, and how the distrubance is calcualted will again be important as both single simulations and co-simulations -are used by ClosedLoopUnitIdentifier to identify the process model including possibly time constants and time-delays. \ No newline at end of file +are used by ClosedLoopUnitIdentifier to identify the process model including possibly time constants and time-delays. + + + + + +### Computational loops other than PID-feedback loops + +The PlantSimulator can deal with computational loops other than PID-feedback loops. These are initalized to steady-state by co-simulating the loop for a number of iterations until the outputs hopefully settle on a steady value. \ No newline at end of file