diff --git a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs index a4d73c3..aa25e9a 100644 --- a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs +++ b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs @@ -268,12 +268,16 @@ public class ClosedLoopUnitIdentifier dataSet.GetTimeBase(), "doDebuggingPlot_disturbanceHLandLF"); dataSet.D = null; - var sim1 = new UnitSimulator(idUnitModelsList[0]); - var sim1results = sim1.Simulate(ref dataSet,false,true); - var sim2 = new UnitSimulator(idUnitModelsList[1]); - var sim2results = sim2.Simulate(ref dataSet, false, true); - var sim3 = new UnitSimulator(idUnitModelsList[2]); - var sim3results = sim3.Simulate(ref dataSet, false, true); + /* var sim1 = new UnitSimulator(idUnitModelsList[0]); + var sim1results = sim1.Simulate(ref dataSet,false,true); + var sim2 = new UnitSimulator(idUnitModelsList[1]); + var sim2results = sim2.Simulate(ref dataSet, false, true); + var sim3 = new UnitSimulator(idUnitModelsList[2]); + var sim3results = sim3.Simulate(ref dataSet, false, true); + */ + (var isOk1,var sim1results) = PlantSimulator.SimulateSingle(dataSet, idUnitModelsList[0], false, 0, false); + (var isOk2, var sim2results) = PlantSimulator.SimulateSingle(dataSet, idUnitModelsList[1], false, 0, false); + (var isOk3, var sim3results) = PlantSimulator.SimulateSingle(dataSet, idUnitModelsList[2], false, 0, false); Plot.FromList( new List { @@ -609,18 +613,21 @@ private Tuple EstimateSISOdisturbanceForProcGain public bool ClosedLoopSim(UnitDataSet unitData, UnitParameters modelParams, PidParameters pidParams, double[] disturbance,string name="") { - // TODO: replace this with a "closed-loop" simulator that uses the PlantSimulator. - if (pidParams == null) { return false; } - var sim = new UnitSimulator(new UnitModel(modelParams)); - unitData.D = disturbance; - var pid = new PidModel(pidParams); + unitData.D = disturbance; + var pidModel = new PidModel(pidParams); + var unitModel = new UnitModel(modelParams); + // TODO: replace this with a "closed-loop" simulator that uses the PlantSimulator. + //var ps = new PlantSimulator(new List { unitModel,pidModel }); + //var inputData = new TimeSeriesDataSet();// TODO: need to map signal names here. + //ps.Simulate(inputData, out var simData); - bool isOk = sim.CoSimulate(pid, ref unitData); + var sim = new UnitSimulator(unitModel); + bool isOk = sim.CoSimulate(pidModel, ref unitData); if (doDebuggingPlot) { Plot.FromList(new List { unitData.Y_sim,unitData.Y_meas, diff --git a/Dynamic/Identification/DisturbanceIdentifier.cs b/Dynamic/Identification/DisturbanceIdentifier.cs index e8100fc..6be78e1 100644 --- a/Dynamic/Identification/DisturbanceIdentifier.cs +++ b/Dynamic/Identification/DisturbanceIdentifier.cs @@ -112,152 +112,6 @@ public class DisturbanceIdentifier { const double numberOfTiConstantsToWaitAfterSetpointChange = 5; - /// - /// Only uses Y_meas and U in unitDataSet, i.e. does not consider feedback - /// - /// - /// - /// - /// - public static DisturbanceIdResult EstDisturbanceBasedOnProcessModel(UnitDataSet unitDataSet, - UnitModel unitModel, int inputIdx = 0) - { - unitModel.WarmStart(); - var sim = new UnitSimulator(unitModel); - unitDataSet.D = null; - double[] y_sim = sim.Simulate(ref unitDataSet); - - DisturbanceIdResult result = new DisturbanceIdResult(unitDataSet); - result.d_est = (new Vec()).Subtract(unitDataSet.Y_meas, y_sim); - - return result; - } - - /// - /// Removes the effect of setpoint and (if relevant any non-pid input) changes from the dataset using the model of pid and unit provided - /// This is the Multiple-input version of this method that is newer and more generl - /// - /// - /// - /// - /// - /// a scrubbed copy of unitDataSet - private static UnitDataSet RemoveSetpointAndOtherInputChangeEffectsFromDataSet_MISO(UnitDataSet unitDataSet, - UnitModel unitModel, int pidInputIdx = 0, PidParameters pidParams = null) - { - if (Vec.IsConstant(unitDataSet.Y_setpoint)) - { - return unitDataSet; - } - - var unitDataSet_setpointAndExternalEffectsRemoved = new UnitDataSet(unitDataSet); - if (unitModel != null && pidParams != null) - { - // co-simulate the unit model output - // with setpoint changes and changes in any external signals, to give a - // "what the output would have been without disturbances" - // that can then be subtracted from the real y_meas to fin and "e" that goes into d_HF - // in disturbance estimation. - - var pidModel1 = new PidModel(pidParams, "PID1"); - var processSim_noDist = new PlantSimulator( - new List { unitModel,pidModel1 }); - processSim_noDist.ConnectModels(unitModel, pidModel1); - processSim_noDist.ConnectModels(pidModel1, unitModel, pidInputIdx); - - var inputData_noDist = new TimeSeriesDataSet(); - var N = unitModel.GetLengthOfInputVector(); - if (unitDataSet.U.GetNColumns() > 1) - { - for (int curColIdx = 0; curColIdx < unitDataSet.U.GetNColumns(); curColIdx++) - { - if (curColIdx == pidInputIdx) - { - continue;// handled by simulating the pid-controller. - } - else - { - var newConstantInput = Vec.Fill(unitDataSet.U.GetColumn(curColIdx)[0], N); - inputData_noDist.Add(processSim_noDist.AddExternalSignal(unitModel, SignalType.External_U, curColIdx), - unitDataSet.U.GetColumn(curColIdx)); - } - } - } - - - inputData_noDist.Add(processSim_noDist.AddExternalSignal(pidModel1, SignalType.Setpoint_Yset),unitDataSet.Y_setpoint ); - inputData_noDist.CreateTimestamps(unitDataSet.GetTimeBase()); - inputData_noDist.SetIndicesToIgnore(unitDataSet.IndicesToIgnore); - var isOk = processSim_noDist.Simulate(inputData_noDist, out TimeSeriesDataSet simData_noDist); - - if (isOk) - { - int idxFirstGoodValue = 0; - if (unitDataSet.IndicesToIgnore != null) - { - if (unitDataSet.GetNumDataPoints() > 0) - { - while (unitDataSet.IndicesToIgnore.Contains(idxFirstGoodValue) && - idxFirstGoodValue < unitDataSet.GetNumDataPoints() - 1) - { - idxFirstGoodValue++; - } - } - } - - var vec = new Vec(); - - // output Y: subtract from Y_meas the Y_meas that results from - var no_disturbance_OutputY = simData_noDist.GetValues(unitModel.GetID(), SignalType.Output_Y); - unitDataSet_setpointAndExternalEffectsRemoved.Y_meas = vec.Subtract(unitDataSet.Y_meas, no_disturbance_OutputY); - - // inputs U - if (unitDataSet_setpointAndExternalEffectsRemoved.U.GetNColumns() > 1) // todo:not general - { - for (int inputIdx = 0; inputIdx < unitDataSet_setpointAndExternalEffectsRemoved.U.GetNColumns(); inputIdx++) - { - var pidOutputU = unitDataSet.U.GetColumn(inputIdx); - var pidDeltaU = vec.Subtract(pidOutputU, pidOutputU[idxFirstGoodValue]); - var newU = vec.Subtract(unitDataSet.U.GetColumn(inputIdx), pidDeltaU); - unitDataSet_setpointAndExternalEffectsRemoved.U = Matrix.ReplaceColumn(unitDataSet_setpointAndExternalEffectsRemoved.U, inputIdx, newU); - } - } - else - { - var pidOutputU = simData_noDist.GetValues(pidModel1.GetID(), SignalType.PID_U); - var pidDeltaU = vec.Subtract(pidOutputU, pidOutputU[idxFirstGoodValue]); - var newU = vec.Subtract(unitDataSet.U.GetColumn(pidInputIdx), pidDeltaU); - unitDataSet_setpointAndExternalEffectsRemoved.U = Matrix.ReplaceColumn(unitDataSet_setpointAndExternalEffectsRemoved.U, pidInputIdx, newU); - } - - // Yset : setpoints(is not used by UnitSimulator anyhow) - unitDataSet_setpointAndExternalEffectsRemoved.Y_setpoint = Vec.Fill(unitDataSet.Y_setpoint[idxFirstGoodValue], unitDataSet.Y_setpoint.Length); - unitDataSet_setpointAndExternalEffectsRemoved.IndicesToIgnore = unitDataSet.IndicesToIgnore; - - /* bool doDebugPlot = false; - if (doDebugPlot) - { - - Shared.EnablePlots(); - Plot.FromList( - new List { - unitDataSet_setpointAndExternalEffectsRemoved.Y_meas, - unitDataSet.Y_meas, - unitDataSet_setpointAndExternalEffectsRemoved.Y_setpoint, - unitDataSet.Y_setpoint, - unitDataSet_setpointAndExternalEffectsRemoved.U.GetColumn(pidInputIdx), - unitDataSet.U.GetColumn(pidInputIdx) - }, - new List { "y1=y_meas(new)", "y1=y_meas(old)", "y1=y_set(new)", "y1=y_set(old)", "y3=u_pid(new)", "y3=u_pid(old)" }, - inputData_noDist.GetTimeBase(), "distIdent_setpointTest"); - Shared.DisablePlots(); - }*/ - } - } - return unitDataSet_setpointAndExternalEffectsRemoved; - } - - /// /// Removes the effect of setpoint and (if relevant any non-pid input) changes from the dataset using the model of pid and unit provided /// @@ -567,9 +421,11 @@ public static DisturbanceIdResult EstimateDisturbance(UnitDataSet unitDataSet_ra // non-disturbance related changes in the dataset producing "unitDataSet_adjusted" var unitDataSet_adjusted = RemoveSetpointAndOtherInputChangeEffectsFromDataSet(unitDataSet_raw, unitModel, pidInputIdx, pidParams); unitModel.WarmStart(); - var sim = new UnitSimulator(unitModel); + // var sim = new UnitSimulator(unitModel); unitDataSet_adjusted.D = null; - double[] y_sim = sim.Simulate(ref unitDataSet_adjusted); + // double[] y_sim = sim.Simulate(ref unitDataSet_adjusted); + (var isOk, var y_sim) = PlantSimulator.SimulateSingle(unitDataSet_adjusted, unitModel, false, 0, false); + if (y_sim == null) { result.zeroReason = DisturbanceSetToZeroReason.UnitSimulatorUnableToRun; diff --git a/Dynamic/Identification/GainSchedIdentifier.cs b/Dynamic/Identification/GainSchedIdentifier.cs index ddc00b4..0616e0f 100644 --- a/Dynamic/Identification/GainSchedIdentifier.cs +++ b/Dynamic/Identification/GainSchedIdentifier.cs @@ -35,9 +35,8 @@ static public class GainSchedIdentifier /// /// /// - static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFittingSpecs gsFittingSpecs = null) + static public GainSchedModel Identify(UnitDataSet dataSet, GainSchedFittingSpecs gsFittingSpecs = null) { - int gainSchedInputIndex = 0; if (gsFittingSpecs != null) { @@ -122,7 +121,7 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting // pass 2: const bool DO_PASS2 = true; const int pass2Width = 0;//0,1 or 2, design parameter about how wide to do pass 2 aroudn pass 1 result.(higher number is at the expense of accuracy) - GainSchedParameters modelToReturn = new GainSchedParameters(); + GainSchedParameters paramsToReturn = new GainSchedParameters(); if (bestModelIdx_pass1 > 1 + pass2Width && bestModelIdx_pass1 < allGainSchedParams.Count() - pass2Width && DO_PASS2) { double? gsSearchMin_pass2 = allGainSchedParams.ElementAt(Math.Max(bestModelIdx_pass1 - 1 - pass2Width, 0)).LinearGainThresholds.First(); @@ -131,14 +130,17 @@ static public GainSchedParameters Identify(UnitDataSet dataSet, GainSchedFitting IdentifyGainScheduledGainsAndSingleThreshold(dataSet, gainSchedInputIndex, true, globalSearchIterationsPass2, gsSearchMin_pass2, gsSearchMax_pass2); (var bestModel_pass2, var bestModelIdx_pass2) = ChooseBestGainScheduledModel(potentialGainschedParametersList_pass2, potentialYsimList_pass2, ref dataSet); - modelToReturn = bestModel_pass2; + paramsToReturn = bestModel_pass2; } else { - modelToReturn = bestModel_pass1; + paramsToReturn = bestModel_pass1; } - EstimateTimeDelay(ref modelToReturn, ref dataSet); - return modelToReturn; + EstimateTimeDelay(ref paramsToReturn, ref dataSet); + paramsToReturn.Fitting = new FittingInfo(); + paramsToReturn.Fitting.WasAbleToIdentify = true; + paramsToReturn.Fitting.SolverID = "Identify(thresholds estimated)"; + return new GainSchedModel(paramsToReturn,"identified"); } /// @@ -219,12 +221,12 @@ private static void EstimateTimeDelay(ref GainSchedParameters gsParams, ref Unit /// the object in which the thresholds to be used are given /// set to false to disable estimation of time delays, default is true /// - public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet, GainSchedFittingSpecs gsFittingSpecs, bool doTimeDelayEstimation = true) + public static GainSchedModel IdentifyForGivenThresholds(UnitDataSet dataSet, GainSchedFittingSpecs gsFittingSpecs, bool doTimeDelayEstimation = true) { var vec = new Vec(); - GainSchedParameters retModel = new GainSchedParameters(); - retModel.GainSchedParameterIndex = gsFittingSpecs.uGainScheduledInputIndex; - retModel.LinearGainThresholds = gsFittingSpecs.uGainThresholds; + GainSchedParameters idParams = new GainSchedParameters(); + idParams.GainSchedParameterIndex = gsFittingSpecs.uGainScheduledInputIndex; + idParams.LinearGainThresholds = gsFittingSpecs.uGainThresholds; // for this to work roubustly, the training set for fitting each model may need to require adding in a "span" of neighboring models. int numberOfInputs = dataSet.U.GetLength(1); double gsVarMinU, gsVarMaxU; @@ -234,7 +236,7 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet var warningNotEnoughExitationBetweenAllThresholds = false; var dataSetCopy = new UnitDataSet(dataSet); // estimate each of the gains one by one - for (int curGainIdx = 0; curGainIdx < retModel.LinearGainThresholds.Count()+1; curGainIdx++) + for (int curGainIdx = 0; curGainIdx < idParams.LinearGainThresholds.Count()+1; curGainIdx++) { double[] uMinFit = new double[numberOfInputs]; double[] uMaxFit = new double[numberOfInputs]; @@ -243,23 +245,23 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet if (curGainIdx == 0) { - gsVarMinU = vec.Min(Array2D.GetColumn(dataSetCopy.U, retModel.GainSchedParameterIndex)); + gsVarMinU = vec.Min(Array2D.GetColumn(dataSetCopy.U, idParams.GainSchedParameterIndex)); } else { - gsVarMinU = retModel.LinearGainThresholds[curGainIdx-1]; + gsVarMinU = idParams.LinearGainThresholds[curGainIdx-1]; } - if (curGainIdx == retModel.LinearGainThresholds.Count()) + if (curGainIdx == idParams.LinearGainThresholds.Count()) { - gsVarMaxU = vec.Max(Array2D.GetColumn(dataSetCopy.U, retModel.GainSchedParameterIndex)); + gsVarMaxU = vec.Max(Array2D.GetColumn(dataSetCopy.U, idParams.GainSchedParameterIndex)); } else { - gsVarMaxU = retModel.LinearGainThresholds[curGainIdx]; + gsVarMaxU = idParams.LinearGainThresholds[curGainIdx]; } for (int idx = 0; idx < numberOfInputs; idx++) { - if (idx == retModel.GainSchedParameterIndex) + if (idx == idParams.GainSchedParameterIndex) { uMinFit[idx] = gsVarMinU; uMaxFit[idx] = gsVarMaxU; @@ -270,10 +272,10 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet uMaxFit[idx] = double.NaN; } } - var idResults = IdentifySingleLinearModelForGivenThresholds(ref dataSetCopy, uMinFit, uMaxFit, retModel.GainSchedParameterIndex,u0,false); + var idResults = IdentifySingleLinearModelForGivenThresholds(ref dataSetCopy, uMinFit, uMaxFit, idParams.GainSchedParameterIndex,u0,false); if (idResults.NotEnoughExitationBetweenAllThresholds) warningNotEnoughExitationBetweenAllThresholds = true; - var uInsideUMinUMax = Vec.GetValuesAtIndices(dataSetCopy.U.GetColumn(retModel.GainSchedParameterIndex), + var uInsideUMinUMax = Vec.GetValuesAtIndices(dataSetCopy.U.GetColumn(idParams.GainSchedParameterIndex), Index.InverseIndices(dataSetCopy.GetNumDataPoints(),dataSetCopy.IndicesToIgnore)); var uMaxObserved = (new Vec()).Max(uInsideUMinUMax); var uMinObserved = (new Vec()).Min(uInsideUMinUMax); @@ -287,34 +289,35 @@ public static GainSchedParameters IdentifyForGivenThresholds(UnitDataSet dataSet } // final gain:above the highest threshold - retModel.LinearGains = linearGains; + idParams.LinearGains = linearGains; // time delay if (gsFittingSpecs.uTimeConstantThresholds == null) { - retModel.TimeConstant_s = new double[] { vec.Mean(timeConstants.ToArray()).Value }; + idParams.TimeConstant_s = new double[] { vec.Mean(timeConstants.ToArray()).Value }; } else if (Vec.Equal(gsFittingSpecs.uTimeConstantThresholds,gsFittingSpecs.uGainThresholds)) { - retModel.TimeConstant_s = timeConstants.ToArray(); - retModel.TimeConstantThresholds = gsFittingSpecs.uTimeConstantThresholds; + idParams.TimeConstant_s = timeConstants.ToArray(); + idParams.TimeConstantThresholds = gsFittingSpecs.uTimeConstantThresholds; } else { - retModel.TimeConstant_s = null;//TODO: currently not supported to find timeconstants separately from gain thresholds. + idParams.TimeConstant_s = null;//TODO: currently not supported to find timeconstants separately from gain thresholds. } - retModel.Fitting = new FittingInfo(); - retModel.Fitting.WasAbleToIdentify = allIdsOk; + idParams.Fitting = new FittingInfo(); + idParams.Fitting.SolverID = "IdentifyForGivenThresholds"; + idParams.Fitting.WasAbleToIdentify = allIdsOk; if (!allIdsOk) - retModel.AddWarning(GainSchedIdentWarnings.UnableToIdentifySomeSubmodels); + idParams.AddWarning(GainSchedIdentWarnings.UnableToIdentifySomeSubmodels); if(warningNotEnoughExitationBetweenAllThresholds) - retModel.AddWarning(GainSchedIdentWarnings.InsufficientExcitationBetweenEachThresholdToBeCertainOfGains); + idParams.AddWarning(GainSchedIdentWarnings.InsufficientExcitationBetweenEachThresholdToBeCertainOfGains); // simulate the model and determine the optimal bias term: - DetermineOperatingPointAndSimulate(ref retModel, ref dataSet); + DetermineOperatingPointAndSimulate(ref idParams, ref dataSet); - EstimateTimeDelay(ref retModel, ref dataSet); + EstimateTimeDelay(ref idParams, ref dataSet); - return retModel; + return new GainSchedModel(idParams,"identified"); } /// diff --git a/Dynamic/Identification/UnitIdentifier.cs b/Dynamic/Identification/UnitIdentifier.cs index 3934941..42b63f4 100644 --- a/Dynamic/Identification/UnitIdentifier.cs +++ b/Dynamic/Identification/UnitIdentifier.cs @@ -118,8 +118,9 @@ public static UnitModel IdentifyLinearDiff(ref UnitDataSet dataSet, FittingSpecs if (model.modelParameters.Fitting.WasAbleToIdentify) { - var simulator = new UnitSimulator(model); - simulator.Simulate(ref dataSet, default, true);// overwrite any y_sim + PlantSimulator.SimulateSingle(dataSet, model, false, 0, true); + //var simulator = new UnitSimulator(model); + //simulator.Simulate(ref dataSet, default, true);// overwrite any y_sim model.SetFittedDataSet(dataSet); } return model; @@ -254,10 +255,9 @@ private static UnitModel Identify_Internal(ref UnitDataSet dataSet, FittingSpecs uNorm = SignificantDigits.Format(uNorm, nDigits); } - TimeSpan span = dataSet.GetTimeSpan(); - double maxExpectedTc_s = span.TotalSeconds / 4; - UnitTimeDelayIdentifier processTimeDelayIdentifyObj = - new UnitTimeDelayIdentifier(dataSet.GetTimeBase(), maxExpectedTc_s); + + + // logic for all curves off an all curves on, treated as special cases bool[] allCurvesDisabled = new bool[u0.Count()]; @@ -275,13 +275,34 @@ private static UnitModel Identify_Internal(ref UnitDataSet dataSet, FittingSpecs var dataset_copy = new UnitDataSet(dataSet); + // find a static model with no time delay or nonlinearity UnitParameters modelParams_StaticAndNoCurvature = EstimateProcessForAGivenTimeDelay (timeDelayIdx, dataSet, false, allCurvesDisabled, FilterTc_s, u0, uNorm, assumeThatYkminusOneApproxXkminusOne); modelList.Add(modelParams_StaticAndNoCurvature); - var warningList = new List(); + // find a dynamic model with no time delay or nonlinearity + // (the time constant gives an upper bound on the time delay) + /* UnitParameters modelParams_DynamicAndNoCurvature = + EstimateProcessForAGivenTimeDelay + (timeDelayIdx, dataSet, true, allCurvesDisabled, + FilterTc_s, u0, uNorm, assumeThatYkminusOneApproxXkminusOne); + modelList.Add(modelParams_DynamicAndNoCurvature); + // modelList.Add(modelParams_DynamicAndNoCurvature); + // double maxExpectedTc_s = Math.Ceiling(modelParams_DynamicAndNoCurvature.TimeConstant_s/ dataSet.GetTimeBase())* dataSet.GetTimeBase() + dataSet.GetTimeBase()*3; + + */ + var warningList = new List(); + + //nb! this is quite a high uppper bound, in the worst case this can cause + // an obscene number of iterations of the below while loop. + TimeSpan span = dataSet.GetTimeSpan(); + double maxExpectedTc_s = span.TotalSeconds / 4; + + UnitTimeDelayIdentifier processTimeDelayIdentifyObj = + new UnitTimeDelayIdentifier(dataSet.GetTimeBase(), maxExpectedTc_s); + ///////////////////////////////////////////////////////////////// // BEGIN WHILE loop to model process for different time delays @@ -422,10 +443,7 @@ private static UnitModel Identify_Internal(ref UnitDataSet dataSet, FittingSpecs { modelParameters.AddWarning(UnitdentWarnings.CorrelatedInputsU); } - modelParameters.FittingSpecs = fittingSpecs; - - // END While loop ///////////////////////////////////////////////////////////////// var model = new UnitModel(modelParameters, dataSet); @@ -433,8 +451,7 @@ private static UnitModel Identify_Internal(ref UnitDataSet dataSet, FittingSpecs // simulate if (modelParameters.Fitting.WasAbleToIdentify) { - var simulator = new UnitSimulator(model); - simulator.Simulate(ref dataSet, default, true);// overwrite any y_sim + PlantSimulator.SimulateSingle(dataSet, model, false, 0, true);// overwrite any y_sim model.SetFittedDataSet(dataSet); } @@ -1124,7 +1141,7 @@ private static void CalculateDynamicUncertainty(RegressionResults regResults,dou /// /// For a given set of paramters and a dataset, find the bias which gives the lowest mean offset - /// (this can be espescially useful when identification uses the "difference" formulation) + /// (this can be especially useful when identification uses the "difference" formulation) /// /// /// @@ -1137,6 +1154,8 @@ static public (double?, double[]) SimulateAndReEstimateBias(UnitDataSet dataSet, var model = new UnitModel(parameters); var simulator = new UnitSimulator(model); var y_sim = simulator.Simulate(ref internalData); + // (var isOk, var y_sim) = PlantSimulator.SimulateSingle(internalData, model,false, 0, true); + var yMeas_exceptIgnoredValues = internalData.Y_meas; var ySim_exceptIgnoredValues = y_sim; if (dataSet.IndicesToIgnore != null) diff --git a/Dynamic/PlantSimulator/PlantSimulator.cs b/Dynamic/PlantSimulator/PlantSimulator.cs index 2af8fd2..500b0ae 100644 --- a/Dynamic/PlantSimulator/PlantSimulator.cs +++ b/Dynamic/PlantSimulator/PlantSimulator.cs @@ -2,6 +2,7 @@ using System.Collections.Generic; using System.Diagnostics; using System.Linq; +using System.Net.Http.Headers; using System.Reflection; using System.Text; using System.Threading; @@ -433,7 +434,7 @@ public bool SimulateSingleInternal(TimeSeriesDataSet inputData, string singleMod } /// - /// Simualte single model to get the output including any additive inputs + /// Simulate single model to get the output including any additive inputs /// /// /// @@ -444,6 +445,72 @@ public bool SimulateSingle(TimeSeriesDataSet inputData, string singleModelName, return SimulateSingle(inputData, singleModelName, false, out simData); } + /// + /// Simulate single model based on a unit data set + /// + /// contains a unit data set that must have U filled, Y_sim will be written here + /// model to simulate + /// if set to true, the simulated result is written to unitData.Y_meas instead of Y_sim + /// if writing to Ymeas, it is possible to add noise of the given amplitude to signal + /// if true, the Y_sim of unitData has the simulation result written two i + /// a tuple, first aa true if able to simulate, otherwise false, second is the simulated time-series + static public (bool, double[]) SimulateSingle(UnitDataSet unitData, UnitModel model,bool writeToYmeas= false, double noiseAmplitude=0, + bool addSimToUnitData=false) + { + var inputData = new TimeSeriesDataSet(); + var singleModelName = "SimulateSingle"; + var modelCopy = new UnitModel(model.GetModelParameters(), 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()); + + PlantSimulator sim = new PlantSimulator(new List { modelCopy }); + + var simData = new TimeSeriesDataSet(); + + var isOk = sim.SimulateSingle(inputData, singleModelName, false, out simData); + + double[] y_sim = simData.GetValues(singleModelName, SignalType.Output_Y); + if (noiseAmplitude > 0) + { + // use a specific seed here, to avoid potential issues with "random unit tests" and not-repeatable + // errors. + Random rand = new Random(1232); + + for (int k = 0; k < y_sim.Count(); k++) + { + y_sim[k] += (rand.NextDouble() - 0.5) * 2 * noiseAmplitude; + } + } + + if (addSimToUnitData) + { + if (writeToYmeas) + { + unitData.Y_meas = y_sim; + } + else + { + unitData.Y_sim = y_sim; + } + } + return (isOk, y_sim); + } + + /// /// Simulate a single model(any ISimulatable model), using inputData as inputs, /// @@ -472,6 +539,8 @@ public bool SimulateSingle(TimeSeriesDataSet inputData, string singleModelName, simData = new TimeSeriesDataSet(); int? N = inputData.GetLength(); + if (N.Value == 0) + return false; int timeIdx = 0; var model = modelDict[singleModelName]; string[] additiveInputIDs = model.GetAdditiveInputIDs(); diff --git a/Dynamic/SimulatableModels/GainSchedModel.cs b/Dynamic/SimulatableModels/GainSchedModel.cs index 263e916..8586ebf 100644 --- a/Dynamic/SimulatableModels/GainSchedModel.cs +++ b/Dynamic/SimulatableModels/GainSchedModel.cs @@ -9,6 +9,7 @@ using TimeSeriesAnalysis.Utility; using System.Reflection; +using Accord.Math; namespace TimeSeriesAnalysis.Dynamic { @@ -64,6 +65,14 @@ public GainSchedModel(GainSchedParameters modelParameters, string ID="not_named" InitSim(modelParameters); } + /// + /// Empty constructor + /// + public GainSchedModel() + { + + } + /// /// Answers if model is simulatable, i.e. has given inputs that are sensible and sufficent. /// @@ -470,7 +479,7 @@ private double CalculateStaticStateWithoutAdditive(double[] inputs, double badVa /// /// optional special value that indicates "Nan" /// - public double? GetSteadyStateOutput(double[] u0, double badDataID) + public double? GetSteadyStateOutput(double[] u0, double badDataID= -9999) { if (modelParameters.LinearGains == null) return 0; @@ -589,10 +598,135 @@ public double[] Iterate(double[] inputs, double timeBase_s,double badValueIndica /// override public string ToString() { - return ""; - } + var writeCulture = new CultureInfo("en-US");// System.Globalization.CultureInfo.InstalledUICulture; + var numberFormat = (System.Globalization.NumberFormatInfo)writeCulture.NumberFormat.Clone(); + numberFormat.NumberDecimalSeparator = "."; + + int sDigits = 3; + int sDigitsUnc = 2; + + int cutOffForUsingDays_s = 86400; + int cutOffForUsingHours_s = 3600; + + StringBuilder sb = new StringBuilder(); + sb.AppendLine(this.GetType().ToString()); + sb.AppendLine("-------------------------"); + if (modelParameters.Fitting == null) + { + sb.AppendLine("a priori model"); + } + else + { + if (modelParameters.Fitting.WasAbleToIdentify) + { + sb.AppendLine("ABLE to identify"); + } + else + { + sb.AppendLine("---NOT able to identify---"); + } + } + + sb.AppendLine("Gain-scheduling parameter index:" + modelParameters.GainSchedParameterIndex); + sb.AppendLine("Operating point U:" + SignificantDigits.Format(modelParameters.OperatingPoint_U, sDigits).ToString(writeCulture) ); + sb.AppendLine("Operating point Y:" + SignificantDigits.Format(modelParameters.OperatingPoint_Y, sDigits).ToString(writeCulture) ); + //////////////////////////////// + // time constant + string timeConstantString = "TimeConstant : "; + + for (int i = 0; i < modelParameters.TimeConstant_s.Count(); i++) + { + if (modelParameters.TimeConstant_s[i] < cutOffForUsingHours_s) + { + timeConstantString += + SignificantDigits.Format(modelParameters.TimeConstant_s[i], sDigits).ToString(writeCulture) + " sec"; + } + else if (modelParameters.TimeConstant_s[i] < cutOffForUsingDays_s) + { + timeConstantString += + SignificantDigits.Format(modelParameters.TimeConstant_s[i] / 3600, sDigits).ToString(writeCulture) + " hours"; + } + else // use days + { + timeConstantString += + SignificantDigits.Format(modelParameters.TimeConstant_s[i] / 86400, sDigits).ToString(writeCulture) + " days"; + } + } + sb.AppendLine(timeConstantString); + //////////////////////////////// + + if (modelParameters.TimeConstantThresholds == null) + sb.AppendLine("Time constant thresholds : [null]"); + else + { + sb.AppendLine("Time constant thresholds : "); + for (int inputIdx = 0; inputIdx < modelParameters.TimeConstantThresholds.Count(); inputIdx++) + { + sb.AppendLine( + "\t" + SignificantDigits.Format(modelParameters.TimeConstantThresholds[inputIdx], sDigits).ToString(writeCulture) + ); + } + } + //////////////////////////////// + // time delay + if (modelParameters.TimeDelay_s < cutOffForUsingHours_s) + { + sb.AppendLine("TimeDelay : " + modelParameters.TimeDelay_s.ToString(writeCulture) + " sec"); + } + else if (modelParameters.TimeDelay_s < cutOffForUsingDays_s) + { + sb.AppendLine("TimeDelay : " + (modelParameters.TimeDelay_s / 3600).ToString(writeCulture) + " sec"); + } + else + { + sb.AppendLine("TimeDelay : " + (modelParameters.TimeDelay_s / 86400).ToString(writeCulture) + " days"); + } + + //////////////////////////////// + sb.AppendLine("Linear gains : "); + for (int inputIdx = 0; inputIdx < modelParameters.LinearGains.Count; inputIdx++) + { + for (int gsVarIdx = 0; gsVarIdx < modelParameters.LinearGains[inputIdx].Count(); gsVarIdx++) + { + sb.AppendLine( + "\t" + SignificantDigits.Format(modelParameters.LinearGains[inputIdx][gsVarIdx], sDigits).ToString(writeCulture) + ); + } + } + + //////////////////////////////// + sb.AppendLine("Linear gains thresholds : "); + for (int inputIdx = 0; inputIdx < modelParameters.LinearGainThresholds.Count(); inputIdx++) + { + sb.AppendLine( + "\t" + SignificantDigits.Format(modelParameters.LinearGainThresholds[inputIdx], sDigits).ToString(writeCulture) + ); + } + //////////////////////////////// + sb.AppendLine("-------------------------"); + if (modelParameters.Fitting != null) + { + /* sb.AppendLine("Fit score(%): " + modelParameters.Fitting.FitScorePrc.ToString(writeCulture)); + + sb.AppendLine("objective(diffs): " + SignificantDigits.Format(modelParameters.Fitting.ObjFunValDiff, 4).ToString(writeCulture)); + sb.AppendLine("R2(diffs): " + SignificantDigits.Format(modelParameters.Fitting.RsqDiff, 4).ToString(writeCulture)); + sb.AppendLine("R2(abs): " + SignificantDigits.Format(modelParameters.Fitting.RsqAbs, 4).ToString(writeCulture)); + */ + sb.AppendLine("model fit data points: " + modelParameters.Fitting.NFittingTotalDataPoints + " of which " + modelParameters.Fitting.NFittingBadDataPoints + " were excluded"); + foreach (var warning in modelParameters.GetWarningList()) + sb.AppendLine("model fit warning :" + warning.ToString()); + if (modelParameters.GetWarningList().Count == 0) + { + sb.AppendLine("model fit : no error or warnings"); + } + + sb.AppendLine("solver: " + modelParameters.Fitting.SolverID); + } + return sb.ToString(); + } + /// /// Warm-starting /// diff --git a/Dynamic/UnitSimulator/UnitSimulator.cs b/Dynamic/UnitSimulator/UnitSimulator.cs index d84e142..9131688 100644 --- a/Dynamic/UnitSimulator/UnitSimulator.cs +++ b/Dynamic/UnitSimulator/UnitSimulator.cs @@ -28,31 +28,6 @@ public UnitSimulator(UnitModel model) this.model = model; } - /// - /// Simulation is written to ymeas instead of ysim. This is useful when creating generic datasets for - /// testing/test driven development. - /// - /// - /// optionally adds noise to the "measured" y (for testing purposes) - public bool SimulateYmeas(ref UnitDataSet processDataSet, double noiseAmplitude=0) - { - if (processDataSet.U == null) - return false; - Simulate(ref processDataSet,true); - if (noiseAmplitude > 0) - { - // use a specific seed here, to avoid potential issues with "random unit tests" and not-repeatable - // errors. - Random rand = new Random(1232); - - for (int k = 0; k < processDataSet.GetNumDataPoints(); k++) - { - processDataSet.Y_meas[k] += (rand.NextDouble()-0.5)*2* noiseAmplitude; - } - } - return true; - } - /// /// Co-simulate a process model and pid-controller(both Y_sim and U_sim) diff --git a/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs b/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs index 2a6603c..299a505 100644 --- a/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs +++ b/TimeSeriesAnalysis.Tests/Examples/SystemIdent.cs @@ -65,11 +65,10 @@ public void NonlinearUnitModel() Bias = bias }; var refModel = new UnitModel(paramtersNoCurvature, "Reference"); - var refSim = new UnitSimulator(refModel); var refData = new UnitDataSet(); refData.U = U; refData.CreateTimeStamps(timeBase_s); - refSim.Simulate(ref refData); + PlantSimulator.SimulateSingle(refData, refModel); // simulate the nonlinear model UnitParameters designParameters = new UnitParameters @@ -83,11 +82,10 @@ public void NonlinearUnitModel() Bias = bias }; var trueModel = new UnitModel(designParameters, "NonlinearModel1"); - var sim = new UnitSimulator(trueModel); var idDataSet = new UnitDataSet(); idDataSet.U = U; idDataSet.CreateTimeStamps(timeBase_s); - sim.SimulateYmeas(ref idDataSet, noiseAmplitude); + PlantSimulator.SimulateSingle(idDataSet, trueModel,true, noiseAmplitude); // do identification of unit model FittingSpecs fittingSpecs = new FittingSpecs(designParameters.U0, designParameters.UNorm); @@ -100,7 +98,7 @@ public void NonlinearUnitModel() new List { "y1=ysim", "y1=ymeas", "y1=yref(linear)", "y3=u1", "y3=u2" }, (int)timeBase_s, "NonlinearUnitModelEx", default); - PlotGain.Plot(idModel, trueModel); + PlotGain.PlotSteadyState(idModel, trueModel); Console.WriteLine(idModel.ToString()); } diff --git a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs b/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs index 7c4b05f..5899412 100644 --- a/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/GainSchedIdentifyTests.cs @@ -25,7 +25,7 @@ private void ConoleOutResult(GainSchedParameters trueParams, GainSchedParameters string delimTxt = " vs "; string accuracy = "F2"; - Console.Write("[estimate]" + delimTxt + "[true]"); + Console.Write("[estimate]" + delimTxt + "[true]"+"\r\n"); for (int i = 0; i < estParams.TimeConstant_s.Length; i++) { @@ -130,11 +130,12 @@ public void FiveGains_CorrectGainsReturned(int ver, int expectedNumWarnings, dou dataSet.Y_meas = simData.GetValues(refModel.ID, SignalType.Output_Y); dataSet.U = Array2D.CreateFromList(new List { inputData.GetValues(refModel.ID,SignalType.External_U)}); dataSet.Times = inputData.GetTimeStamps(); - var gsParams = GainSchedIdentifier.IdentifyForGivenThresholds(dataSet, gsFittingSpecs); + var gsModel = GainSchedIdentifier.IdentifyForGivenThresholds(dataSet, gsFittingSpecs); // console out + ConoleOutResult(refParams, gsModel.GetModelParameters()); - ConoleOutResult(refParams, gsParams); + Console.WriteLine(gsModel); // Plotting gains(debugging) bool doPlots = false; @@ -149,23 +150,22 @@ public void FiveGains_CorrectGainsReturned(int ver, int expectedNumWarnings, dou timeBase_s, "GainEstOnly_CorrectGainsReturned"); - GainSchedModel gsModel = new GainSchedModel(gsParams,"ident_model"); gsModel.SetOutputID("y_meas"); gsModel.SetInputIDs((new List { "u_1" }).ToArray()); GainSchedModel referenceModel = new GainSchedModel(refParams,"ref_model"); referenceModel.SetInputIDs(gsModel.GetModelInputIDs()); referenceModel.SetOutputID(gsModel.GetOutputID()); - - PlotGain.Plot(gsModel, referenceModel); + PlotGain.PlotSteadyState(gsModel, referenceModel, "steady state gains", new double[] { 0}, new double[] { 15}); + PlotGain.PlotGainSched(gsModel, referenceModel,"gain-scheduling"); Shared.DisablePlots(); } // asserts - Assert.IsTrue(gsParams.Fitting.WasAbleToIdentify); - Assert.AreEqual(expectedNumWarnings, gsParams.GetWarningList().Count()); + Assert.IsTrue(gsModel.GetModelParameters().Fitting.WasAbleToIdentify); + Assert.AreEqual(expectedNumWarnings, gsModel.GetModelParameters().GetWarningList().Count()); for (int i = 0; i < refParams.LinearGains.Count; i++) { - DiffLessThan(refParams.LinearGains[i][0], gsParams.LinearGains[i][0], gainTolerancePrc,i); + DiffLessThan(refParams.LinearGains[i][0], gsModel.GetModelParameters().LinearGains[i][0], gainTolerancePrc,i); } } @@ -227,7 +227,7 @@ public void TimeDelay_TDEstOk(int timeDelaySamples) unitData.Y_meas = (new Vec()).Add(Vec.Rand(simY1.Length, -noiseAmp, noiseAmp, 454), simY1); // Act - GainSchedParameters idParams = GainSchedIdentifier.Identify(unitData); + var idModel = GainSchedIdentifier.Identify(unitData); // plot bool doPlot = false; @@ -244,8 +244,8 @@ public void TimeDelay_TDEstOk(int timeDelaySamples) Shared.DisablePlots(); } - ConoleOutResult(trueGSparams, idParams); - Assert.That(Math.Abs(idParams.TimeDelay_s - trueGSparams.TimeDelay_s) { est_model }); - inputData.Add(estPlantSim.AddExternalSignal(est_model, SignalType.External_U, (int)INDEX.FIRST), u); + var estPlantSim = new PlantSimulator(new List { idModel }); + inputData.Add(estPlantSim.AddExternalSignal(idModel, SignalType.External_U, (int)INDEX.FIRST), u); var IdentifiedisSimulatable = estPlantSim.Simulate(inputData, out TimeSeriesDataSet IdentifiedsimData); - double[] simY2 = IdentifiedsimData.GetValues(est_model.GetID(), SignalType.Output_Y); + double[] simY2 = IdentifiedsimData.GetValues(idModel.GetID(), SignalType.Output_Y); // plot bool doPlot = false; if (doPlot) @@ -479,21 +478,21 @@ public void TwoGains_TCAndThresholdFoundOk(double TimeConstant1_s, Shared.DisablePlots(); } - ConoleOutResult(trueParams, est_params); + ConoleOutResult(trueParams, idModel.GetModelParameters()); // Asserts - Assert.IsTrue(Math.Abs(est_params.LinearGainThresholds.First() - trueParams.LinearGainThresholds.First()) < threshold_tol); + Assert.IsTrue(Math.Abs(idModel.GetModelParameters().LinearGainThresholds.First() - trueParams.LinearGainThresholds.First()) < threshold_tol); SISOTests.CommonAsserts(inputData, IdentifiedsimData, estPlantSim); - int min_number_of_time_constants = Math.Min(est_params.TimeConstant_s.Length, trueParams.TimeConstant_s.Length); + int min_number_of_time_constants = Math.Min(idModel.GetModelParameters().TimeConstant_s.Length, trueParams.TimeConstant_s.Length); for (int k = 0; k < min_number_of_time_constants; k++) { - Assert.That(est_params.TimeConstant_s[k].IsGreaterThanOrEqual(Math.Max(0,Math.Min(TimeConstant1_s,TimeConstant2_s) - TimeConstantAllowedDev_s)), + Assert.That(idModel.GetModelParameters().TimeConstant_s[k].IsGreaterThanOrEqual(Math.Max(0,Math.Min(TimeConstant1_s,TimeConstant2_s) - TimeConstantAllowedDev_s)), "Too low time constant " + k.ToString()); - Assert.That(est_params.TimeConstant_s[k].IsLessThanOrEqual(Math.Max(TimeConstant1_s, TimeConstant2_s) + TimeConstantAllowedDev_s), + Assert.That(idModel.GetModelParameters().TimeConstant_s[k].IsLessThanOrEqual(Math.Max(TimeConstant1_s, TimeConstant2_s) + TimeConstantAllowedDev_s), "Too high time constant " + k.ToString()); } - Assert.That(Math.Abs(est_params.OperatingPoint_Y - trueParams.OperatingPoint_Y) < operatingy_tol); + Assert.That(Math.Abs(idModel.GetModelParameters().OperatingPoint_Y - trueParams.OperatingPoint_Y) < operatingy_tol); } diff --git a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs b/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs index db8715d..7840531 100644 --- a/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/PlantSimulatorSISOTests.cs @@ -138,14 +138,31 @@ public void SimulateSingle_NullGains_RunsWithZeroOutput() // now test that "simulateSingle" produces the same result! var isOk2 = plantSim.SimulateSingle(inputData, processModel1.ID, out TimeSeriesDataSet simData2); - /* - Plot.FromList(new List { + + + //plots + + bool doPlot = false; + if (doPlot) + { + Shared.EnablePlots(); + Plot.FromList(new List { simData.GetValues(processModel1.GetID(),SignalType.Output_Y), simData2.GetValues(processModel1.GetID(),SignalType.Output_Y), inputData.GetValues(processModel1.GetID(),SignalType.External_U)}, - new List { "y1=y_sim1", "y1=y_sim1(v2)", "y3=u" }, - timeBase_s, "UnitTest_SingleSISO"); - */ + new List { "y1=y_sim1", "y1=y_sim1(v2)", "y3=u" }, + timeBase_s, "UnitTest_SingleSISO"); + Shared.DisablePlots(); + } + + + // asserts + + Assert.IsTrue(isOk2); + + Assert.AreEqual(simData.GetValues(processModel1.GetID(), SignalType.Output_Y), + simData2.GetValues(processModel1.GetID(), SignalType.Output_Y)); + } diff --git a/TimeSeriesAnalysis.Tests/Tests/SysIDUnitTests.cs b/TimeSeriesAnalysis.Tests/Tests/SysIDUnitTests.cs index 2b4143d..ff7bb95 100644 --- a/TimeSeriesAnalysis.Tests/Tests/SysIDUnitTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/SysIDUnitTests.cs @@ -36,13 +36,21 @@ public void UnitSimulate_TimeDelay(int timeDelay_s) UnitDataSet dataSet = new UnitDataSet(); dataSet.U = Array2D.CreateFromList(new List { u1 }); dataSet.CreateTimeStamps(timeBase_s); - var simulator = new UnitSimulator(model); - var ret = simulator.Simulate(ref dataSet); - // Plot.FromList(new List{ dataSet.Y_sim,u1},new List{"y1=ymeas ","y3=u1"}, timeBase_s); - Assert.IsNotNull(ret); - Assert.IsTrue(dataSet.Y_sim[30+ timeDelay_s] == 0,"step should not arrive at y_sim too early"); - Assert.IsTrue(dataSet.Y_sim[31+ timeDelay_s] == 1, "steps should be delayed exactly timeDelay_s later "); + (var isOk,var 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 "); } } @@ -70,10 +78,11 @@ public UnitDataSet CreateDataSet(UnitParameters designParameters, double[,] U, d dataSet.U = U; dataSet.BadDataID = badValueId; dataSet.CreateTimeStamps(timeBase_s); - var simulator = new UnitSimulator(model); + // var simulator = new UnitSimulator(model); if (doNonWhiteNoise) { - simulator.SimulateYmeas(ref dataSet, 0); + PlantSimulator.SimulateSingle(dataSet, model, true, 0,true); + // simulator.SimulateYmeas(ref dataSet, 0); double rand = 0; var randObj = new Random(45466545); for (int i = 0; i < dataSet.Y_meas.Length; i++) @@ -84,7 +93,8 @@ public UnitDataSet CreateDataSet(UnitParameters designParameters, double[,] U, d } else { - simulator.SimulateYmeas(ref dataSet, noiseAmplitude); + PlantSimulator.SimulateSingle(dataSet, model, true, noiseAmplitude,true); + // simulator.SimulateYmeas(ref dataSet, noiseAmplitude); } if (addInBadDataToYmeasAndU) @@ -634,6 +644,68 @@ public void YminFit_IsExcludedOk(double bias, double timeConstant_s, int timeDel DefaultAsserts(model, designParameters); } + /* + [TestCase(0.10, 5, 5, Category = "Nonlinear,Delayed,Dynamic")] + + + public void I1_NonLinear_RampDown(double curvature, double timeConstant_s, int timeDelay_s) + { + double bias = 1; + double noiseAmplitude = 0.01; + + double[] u1 = TimeSeriesCreator.Ramp(60,80,20,20,20); + double[,] U = Array2D.CreateFromList(new List { u1 }); + + UnitParameters designParameters = new UnitParameters + { + TimeConstant_s = timeConstant_s, + TimeDelay_s = timeDelay_s, + LinearGains = new double[] { 3 }, + Curvatures = new double[] { curvature }, + UNorm = new double[] { 1.1 }, + U0 = new double[] { 50 }, + Bias = bias + }; + + UnitParameters paramtersNoCurvature = new UnitParameters + { + TimeConstant_s = timeConstant_s, + TimeDelay_s = timeDelay_s, + LinearGains = designParameters.LinearGains, + UNorm = designParameters.UNorm, + U0 = designParameters.U0, + Bias = bias + }; + var fittingSpecs = new FittingSpecs(designParameters.U0, designParameters.UNorm); + var refModel = new UnitModel(paramtersNoCurvature, "reference"); + + var sim = new PlantSimulator(new List { refModel }); + var inputData = new TimeSeriesDataSet(); + inputData.Add(sim.AddExternalSignal(refModel, SignalType.External_U), u1); + inputData.CreateTimestamps(timeBase_s); + var isOk = sim.Simulate(inputData, out TimeSeriesDataSet refData); + + var model = CreateDataAndIdentify(designParameters, U, timeBase_s, fittingSpecs, noiseAmplitude); + model.ID = "fitted"; + string caseId = TestContext.CurrentContext.Test.Name; + + //TODO: disable plotting,when done1!!! + Shared.EnablePlots(); + Plot.FromList(new List { model.GetFittedDataSet().Y_sim, + model.GetFittedDataSet().Y_meas, + refData.GetValues(refModel.GetID(),SignalType.Output_Y), + u1 }, + new List { "y1=ysim", "y1=ymeas", "y1=yref(linear)", "y3=u1" }, (int)timeBase_s, caseId, default, + caseId.Replace("(", "").Replace(")", "").Replace(",", "_")); + + PlotGain.Plot(model, new UnitModel(designParameters,"Design"), "RampDownUnitTest"); + + Shared.DisablePlots(); + + DefaultAsserts(model, designParameters); + } + */ + [TestCase(0.4, 0, 0, Category = "Nonlinear")] [TestCase(0.2, 0, 0, Category = "Nonlinear")] @@ -656,7 +728,7 @@ public void YminFit_IsExcludedOk(double bias, double timeConstant_s, int timeDel [TestCase(-0.2, 5, 5, Category = "Nonlinear,Delayed,Dynamic")] [TestCase(-0.4, 5, 5, Category = "Nonlinear,Delayed,Dynamic")] - public void I1_NonLinear(double curvature, double timeConstant_s, int timeDelay_s) + public void I1_NonLinear_Threesteps(double curvature, double timeConstant_s, int timeDelay_s) { double bias = 1; double noiseAmplitude = 0.01; @@ -710,7 +782,7 @@ public void I1_NonLinear(double curvature, double timeConstant_s, int timeDelay_ [TestCase(-0.4, 5, 0, Category = "Nonlinear,Dynamic")] [TestCase(0.4, 5, 5, Category = "Nonlinear,Delayed,Dynamic")] [TestCase(0.2, 5, 5, Category = "Nonlinear,Delayed,Dynamic")] - public void I2_OneNonlinearInput(double curvature, double timeConstant_s, int timeDelay_s) + public void I2_OneNonlinearInput_SixSteps(double curvature, double timeConstant_s, int timeDelay_s) { var model = I2_Internal(curvature,timeConstant_s,timeDelay_s,false); @@ -727,7 +799,7 @@ public void I2_OneNonlinearInput(double curvature, double timeConstant_s, int ti [TestCase(0.4, 5, 5, Category = "Nonlinear,Delayed,Dynamic")] [TestCase(-0.4, 5, 5, Category = "Nonlinear,Delayed,Dynamic")] - public void I2_NonLinear(double curvature, double timeConstant_s, int timeDelay_s,bool curvatureOnBothInputs=true) + public void I2_NonLinear_SixSteps(double curvature, double timeConstant_s, int timeDelay_s,bool curvatureOnBothInputs=true) { I2_Internal(curvature, timeConstant_s, timeDelay_s, curvatureOnBothInputs); } @@ -801,7 +873,7 @@ UnitModel I2_Internal(double curvature, double timeConstant_s, int timeDelay_s, [TestCase(1,15,0, Category ="Dynamic")] [TestCase(1, 0, 2, Category = "Delayed")] - public void I2_Linear(double bias, double timeConstant_s, int timeDelay_s) + public void I2_Linear_Twosteps(double bias, double timeConstant_s, int timeDelay_s) { double noiseAmplitude = 0.01; double[] u1 = TimeSeriesCreator.Step(50, 100, 0, 1); @@ -923,7 +995,7 @@ public void ExcludeDisturbance_I2_Linear(double bias, double timeConstant_s, int [TestCase(1, 0, Category = "Static")] [TestCase(0, 20, Category = "Dynamic")] [TestCase(1, 20, Category = "Dynamic")] - public void I3_Linear_singlesteps(double bias, double timeConstant_s) + public void I3_Linear_OneStep(double bias, double timeConstant_s) { double noiseAmplitude = 0.01; double[] u1 = TimeSeriesCreator.Step(50, 100 ,0,1) ; diff --git a/TimeSeriesAnalysis.csproj b/TimeSeriesAnalysis.csproj index 905476c..f5225ce 100644 --- a/TimeSeriesAnalysis.csproj +++ b/TimeSeriesAnalysis.csproj @@ -14,7 +14,7 @@ False https://github.com/equinor/TimeSeriesAnalysis.git readme.md - 1.3.18 + 1.3.19 Equinor Equinor true diff --git a/Utility/PlotGain.cs b/Utility/PlotGain.cs index fac656c..fa52994 100644 --- a/Utility/PlotGain.cs +++ b/Utility/PlotGain.cs @@ -16,19 +16,23 @@ public class PlotGain { const int numberOfPlotPoints = 30; + private static string ReplaceCommentIllegalStrings(string comment) + { + if (comment == null) + return null; + return comment.Replace(" ","_").Replace(";",""); + } + /// /// Plots an "x-y" plot of the steady-state gains of one or two models. /// /// model of gains to be plotted /// optional seond model to be compared in the plots + /// comment to be added to figure /// optional umin array over which to plot gain plots /// optional umax aray over which to plot gain plots - public static void Plot(UnitModel model1, UnitModel model2 = null, double[] uMin=null, double[] uMax=null) + public static void PlotSteadyState(UnitModel model1, UnitModel model2 = null, string comment = null,double[] uMin=null, double[] uMax=null) { - // TODO: could also plot the uminfit/umaxfit if applicable - //double[] uMinFit = model.modelParameters.FittingSpecs.U_min_fit; - //double[] uMaxFit = model.modelParameters.FittingSpecs.U_max_fit; - string outputId = model1.outputID; if (outputId == null) { @@ -53,14 +57,14 @@ public static void Plot(UnitModel model1, UnitModel model2 = null, double[] uMin List tables_m1 = new List(); if (uMin != null && uMax != null) { - tables_m1 = CreateXYTablesFromUnitModel(model1, model1Name, inputIdx, uMin, uMax); + tables_m1 = CreateSteadyStateXYTablesFromUnitModel(model1, model1Name, inputIdx, uMin, uMax); } else if (model1.modelParameters.Fitting == null) { if (model2 != null) { if (model2.modelParameters.Fitting != null) - tables_m1 = CreateXYTablesFromUnitModel(model1, model1Name, inputIdx, + tables_m1 = CreateSteadyStateXYTablesFromUnitModel(model1, model1Name, inputIdx, model2.modelParameters.Fitting.Umin, model2.modelParameters.Fitting.Umax); else throw new Exception("Unable to plot, no FittingInfo in model(s),specify Umin/Umax and rerun."); @@ -68,7 +72,7 @@ public static void Plot(UnitModel model1, UnitModel model2 = null, double[] uMin } else { - tables_m1 = CreateXYTablesFromUnitModel(model1, model1Name, inputIdx); + tables_m1 = CreateSteadyStateXYTablesFromUnitModel(model1, model1Name, inputIdx); } if (model2 == null) @@ -85,39 +89,42 @@ public static void Plot(UnitModel model1, UnitModel model2 = null, double[] uMin List tables_m2 = new List(); if (uMin != null && uMax != null) { - tables_m2 = CreateXYTablesFromUnitModel(model2, model2Name, inputIdx, uMin, uMax); + tables_m2 = CreateSteadyStateXYTablesFromUnitModel(model2, model2Name, inputIdx, uMin, uMax); } else if (model2.modelParameters.Fitting == null) { if (model1.modelParameters.Fitting != null) - tables_m2 = CreateXYTablesFromUnitModel(model2, model1Name, inputIdx, + tables_m2 = CreateSteadyStateXYTablesFromUnitModel(model2, model1Name, inputIdx, model1.modelParameters.Fitting.Umin, model1.modelParameters.Fitting.Umax); else throw new Exception("Unable to plot, no FittingInfo in model(s),specify Umin/Umax and rerun."); } else - tables_m2 = CreateXYTablesFromUnitModel(model2, model2Name, inputIdx); + tables_m2 = CreateSteadyStateXYTablesFromUnitModel(model2, model2Name, inputIdx); var plotTables = new List(); plotTables.AddRange(tables_m1); plotTables.AddRange(tables_m2); - PlotXY.FromTables(plotTables, outputId + "_" + inputSignalID + "_gains"); + PlotXY.FromTables(plotTables, outputId + "_" + inputSignalID + "_gains", comment); } } + /// - /// Plots an "x-y" plot of the steady-state gains of one or two models. + /// Plots an "x-y" plot of the steady-state of one or two gain-scheduled models. /// /// model of gains to be plotted /// optional seond model to be compared in the plots /// comment to be added to figure + /// optional umin array over which to plot gain plots + /// optional umax array over which to plot gain plots - public static void Plot(GainSchedModel model1, GainSchedModel model2 = null, string comment= null) + public static void PlotSteadyState(GainSchedModel model1, GainSchedModel model2 = null, string comment = null, + double[] uMin = null, double[] uMax = null) { - // TODO: could also plot the uminfit/umaxfit if applicable - //double[] uMinFit = model.modelParameters.FittingSpecs.U_min_fit; - //double[] uMaxFit = model.modelParameters.FittingSpecs.U_max_fit; + comment = ReplaceCommentIllegalStrings(comment); + string outputId = model1.outputID; if (outputId == null) @@ -142,38 +149,83 @@ public static void Plot(GainSchedModel model1, GainSchedModel model2 = null, str } List tables_m1 = new List(); - double[] uMin = null; - double[] uMax = null; - if (model1.modelParameters.Fitting != null) + if (model1.modelParameters.Fitting != null && uMin == null && uMax == null) { - uMin = model1.modelParameters.Fitting.Umin; - uMax = model1.modelParameters.Fitting.Umax; + uMin = model1.modelParameters.Fitting.Umin; + uMax = model1.modelParameters.Fitting.Umax; } - tables_m1 = CreateXYTablesFromGainSchedModel(model1, model1Name, inputIdx, uMin, uMax); + tables_m1 = CreateSteadyStateXYTablesFromGainSchedModel(model1, model1Name, inputIdx,uMin,uMax); - /*if (uMin != null && uMax != null) + if (model2 == null) { - tables_m1 = CreateXYTablesFromModel(model1, model1Name, inputIdx, uMin, uMax); + PlotXY.FromTables(tables_m1, "ss" +outputId + "_" + inputSignalID + "_gains", comment); + return; } - else if (model1.modelParameters.Fitting == null) + ////////////////////////////////// + // if model2 is not null, then plot it as well + if (model2.ID != null && model2.ID != "not_named") { - if (model2 != null) - { - if (model2.modelParameters.Fitting != null) - tables_m1 = CreateXYTablesFromModel(model1, model1Name, inputIdx, - model2.modelParameters.Fitting.Umin, model2.modelParameters.Fitting.Umax); - else - throw new Exception("Unable to plot, no FittingInfo in model(s),specify Umin/Umax and rerun."); - } + model2Name = model2.ID; } - else + List tables_m2 = CreateSteadyStateXYTablesFromGainSchedModel(model2, model2Name, inputIdx,uMin,uMax); + + var plotTables = new List(); + plotTables.AddRange(tables_m1); + plotTables.AddRange(tables_m2); + PlotXY.FromTables(plotTables, outputId + "_" + inputSignalID + "_gains", comment); + } + } + + + + /// + /// Plots an "x-y" plot of the gain-scheduled "linear gains" of one or two gain-scheduled models. + /// + /// Note that it makes no sense to give the "min or "max" of these plots + /// + /// + /// model of gains to be plotted + /// optional seond model to be compared in the plots + /// comment to be added to figure + + + public static void PlotGainSched(GainSchedModel model1, GainSchedModel model2 = null, string comment= null) + { + comment = ReplaceCommentIllegalStrings(comment); + + string outputId = model1.outputID; + if (outputId == null) + { + outputId = "output"; + } + + var model1Name = "model1"; + var model2Name = "model2"; + + if (model1.ID != null && model1.ID != "not_named") + { + model1Name = model1.ID; + } + + for (var inputIdx = 0; inputIdx < model1.ModelInputIDs.Count(); inputIdx++) + { + var inputSignalID = model1.ModelInputIDs[inputIdx]; + if (inputSignalID == null) { - tables_m1 = CreateXYTablesFromModel(model1, model1Name, inputIdx); + inputSignalID = "input" + inputIdx; + } + List tables_m1 = new List(); + + /* if (model1.modelParameters.Fitting != null && uMin == null && uMax == null) + { + uMin = model1.modelParameters.Fitting.Umin; + uMax = model1.modelParameters.Fitting.Umax; }*/ + tables_m1 = CreateGainXYTablesFromGainSchedModel(model1, model1Name, inputIdx); if (model2 == null) { - PlotXY.FromTables(tables_m1, outputId + "_" + inputSignalID + "_gains", comment); + PlotXY.FromTables(tables_m1, "gs_" + outputId + "_" + inputSignalID + "_gains", comment); return; } ////////////////////////////////// @@ -182,23 +234,7 @@ public static void Plot(GainSchedModel model1, GainSchedModel model2 = null, str { model2Name = model2.ID; } - /*List tables_m2 = new List(); - if (uMin != null && uMax != null) - { - tables_m2 = CreateXYTablesFromGainSchedModel(model2, model2Name, inputIdx, uMin, uMax); - } - else if (model2.modelParameters.Fitting == null) - { - if (model1.modelParameters.Fitting != null) - tables_m2 = CreateXYTablesFromGainSchedModel(model2, model1Name, inputIdx, - model1.modelParameters.Fitting.Umin, model1.modelParameters.Fitting.Umax); - else - throw new Exception("Unable to plot, no FittingInfo in model(s),specify Umin/Umax and rerun."); - - } - else - tables_m2 = CreateXYTablesFromGainSchedModel(model2, model2Name, inputIdx);*/ - List tables_m2 = CreateXYTablesFromGainSchedModel(model2, model2Name, inputIdx, uMin, uMax); + List tables_m2 = CreateGainXYTablesFromGainSchedModel(model2, model2Name, inputIdx); var plotTables = new List(); plotTables.AddRange(tables_m1); @@ -207,11 +243,24 @@ public static void Plot(GainSchedModel model1, GainSchedModel model2 = null, str } } - - - static private List CreateXYTablesFromGainSchedModel(GainSchedModel model, string modelName, int inputIdx, + /// + /// This method creates a table of the steady-state calculated model + /// + /// + /// + /// + /// + /// + /// + static private List CreateSteadyStateXYTablesFromGainSchedModel(GainSchedModel model, string modelName, int inputIdx, double[] uMinExt = null, double[] uMaxExt = null) { + + if (inputIdx > 0) + { + Console.WriteLine("currently only supports gain-sched variables with one input"); + return null; + } string outputId = model.outputID; if (outputId == null) { @@ -229,78 +278,35 @@ static private List CreateXYTablesFromGainSchedModel(GainSchedModel mod inputSignalID = "input" + inputIdx; } - XYTable xyTableGain = new XYTable("gain_" + modelName + "_fit", new List { inputSignalID, outputId }, XYlineType.line); - //XYTable xyTableU0 = new XYTable("u0_" + modelName, new List { inputSignalID, outputId }, XYlineType.withMarkers); + XYTable xyTableGain = new XYTable("ssgain_" + modelName + "_fit", new List { inputSignalID, outputId }, XYlineType.line); + XYTable xyTableU0 = new XYTable("op_" + modelName, new List { inputSignalID, outputId }, XYlineType.withMarkers); + xyTableU0.AddRow(new double[] { model.GetModelParameters().OperatingPoint_U, model.GetModelParameters().OperatingPoint_Y }); - var thresholdIdx = -1; - for (int i = 0; i < model.modelParameters.LinearGains.Count(); i++ ) + double uMin = 0, uMax = 100; + if (uMinExt != null) + uMin = uMinExt[inputIdx]; + if (uMaxExt != null) + uMax = uMaxExt[inputIdx]; + + // TODO: this is not completely general, if gain-sched model has more than one input + var u0 = model.modelParameters.OperatingPoint_U; + // double[] u_vec = new double[u0.Length]; + double[] u_vec = new double[1]; +// Array.Copy(u0, u_vec, u0.Length); + for (double uCurInputCurVal = uMin; uCurInputCurVal < uMax; uCurInputCurVal += (uMax - uMin) / numberOfPlotPoints) { - var gain = model.modelParameters.LinearGains.ElementAt(i); - // note that usually there is one less threshold than there are gains, - - double x = 0; - if (thresholdIdx == -1) - { - if (uMinExt == null) - { - x = 0;// todo: this is not completely general - } - else - { - x = uMinExt.ElementAt(inputIdx); - } - } - else - x = model.modelParameters.LinearGainThresholds.ElementAt(thresholdIdx); - xyTableGain.AddRow(new double[] { x,gain.ElementAt(inputIdx) } ); - //xyTableGain.AddRow(new double[] { gain.ElementAt(inputIdx),x }); - thresholdIdx++; + u_vec[inputIdx] = uCurInputCurVal; + var y = model.GetSteadyStateOutput(u_vec); + if (y.HasValue) + xyTableGain.AddRow(new double[] { uCurInputCurVal, y.Value }); } - - /* { - var fittingInfo = model.modelParameters.Fitting; - double uMax = 100, uMin = 0; - if (fittingInfo == null) - { - if (uMinExt != null && uMaxExt != null) - { - uMin = uMinExt.ElementAt(inputIdx); - uMax = uMaxExt.ElementAt(inputIdx); - } - else - { - throw new Exception("no Umin/Umax given externall or found in fittingInfo."); - } - } - - else - { - uMax = fittingInfo.Umax.ElementAt(inputIdx); - uMin = fittingInfo.Umin.ElementAt(inputIdx); - } - var u0 = model.modelParameters.U0; - var y0 = model.GetSteadyStateOutput(u0); - if (y0.HasValue) - xyTableU0.AddRow(new double[] { u0.ElementAt(inputIdx), y0.Value }); - - double[] u_vec = new double[u0.Length]; - Array.Copy(u0, u_vec, u0.Length); - for (double uCurInputCurVal = uMin; uCurInputCurVal < uMax; uCurInputCurVal += (uMax - uMin) / numberOfPlotPoints) - { - u_vec[inputIdx] = uCurInputCurVal; - var y = model.GetSteadyStateOutput(u_vec); - if (y.HasValue) - xyTableGain.AddRow(new double[] { uCurInputCurVal, y.Value }); - } - }*/ - return new List { xyTableGain }; + return new List { xyTableGain,xyTableU0}; } + - - - static private List CreateXYTablesFromUnitModel(UnitModel model, string modelName, int inputIdx, + static private List CreateSteadyStateXYTablesFromUnitModel(UnitModel model, string modelName, int inputIdx, double[] uMinExt=null, double[] uMaxExt=null) { string outputId = model.outputID; @@ -320,7 +326,7 @@ static private List CreateXYTablesFromUnitModel(UnitModel model, string inputSignalID = "input" + inputIdx; } - XYTable xyTableGain = new XYTable("gain_" + modelName + "_fit", new List { inputSignalID, outputId }, XYlineType.line); + XYTable xyTableGain = new XYTable("ssgain_" + modelName + "_fit", new List { inputSignalID, outputId }, XYlineType.line); XYTable xyTableU0= new XYTable("u0_" + modelName, new List { inputSignalID, outputId }, XYlineType.withMarkers); { var fittingInfo = model.modelParameters.Fitting; @@ -361,7 +367,65 @@ static private List CreateXYTablesFromUnitModel(UnitModel model, string return new List { xyTableGain, xyTableU0 }; } + /////////////////////////////////////// + /// + //// "gain" table + + /// + /// This method creates a table from the LinearGains contained in model + /// + /// + /// + /// + /// + static private List CreateGainXYTablesFromGainSchedModel(GainSchedModel model, string modelName, int inputIdx) + { + string outputId = model.outputID; + if (outputId == null) + { + outputId = "output"; + } + + if (model.ID != null && model.ID != "not_named") + { + modelName = model.ID; + } + + var inputSignalID = model.ModelInputIDs[inputIdx]; + if (inputSignalID == null) + { + inputSignalID = "input" + inputIdx; + } + + XYTable xyTableGain = new XYTable("gsgain_" + modelName + "_fit", new List { inputSignalID, outputId }, XYlineType.line); + + var thresholdIdx = -1; + + for (int i = 0; i < model.modelParameters.LinearGains.Count(); i++) + { + var gain = model.modelParameters.LinearGains.ElementAt(i); + // note that usually there is one less threshold than there are gains, + double x = 0; + if (thresholdIdx == -1) + { + // if (uMinExt == null) + { + x = 0;// todo: this is not completely general + } + /* else + { + x = uMinExt.ElementAt(inputIdx); + }*/ + } + else + x = model.modelParameters.LinearGainThresholds.ElementAt(thresholdIdx); + xyTableGain.AddRow(new double[] { x, gain.ElementAt(inputIdx) }); + thresholdIdx++; + } + + return new List { xyTableGain }; + }