diff --git a/Dynamic/Identification/ClosedLoopGainGlobalSearchResults.cs b/Dynamic/Identification/ClosedLoopGainGlobalSearchResults.cs index 657fb82..e18919c 100644 --- a/Dynamic/Identification/ClosedLoopGainGlobalSearchResults.cs +++ b/Dynamic/Identification/ClosedLoopGainGlobalSearchResults.cs @@ -1,4 +1,5 @@ using Accord.Math; +using Accord.Statistics; using System; using System.Collections.Generic; using System.Globalization; @@ -56,20 +57,56 @@ internal ClosedLoopGainGlobalSearchResults() uPidDistanceList = new List(); dEstDistanceList = new List(); covBtwDestAndYsetList = new List(); - pidLinearProcessGainList = new List(); covBtwDestAndUexternal = new List(); } - public void Add(double gain, UnitParameters unitParameters, double covBtwDestAndYset, - double uPidDistance,double dEstDistance, - double covBtwDestAndUexternal) + /// + /// Add the result of one paramter study as part of the global search + /// + /// paramters that have been studied + /// + /// + /// + /// + internal void Add(UnitParameters unitParameters, UnitDataSet dataSet, double[] dEst, double[] u_pid_adjusted, int pidInputIdx) { - pidLinearProcessGainList.Add(gain); + var vec = new Vec(dataSet.BadDataID); + + double covarianceBtwDestAndYset = Math.Abs(CorrelationCalculator.CorrelateTwoVectors( + dEst, dataSet.Y_setpoint, dataSet.IndicesToIgnore)); + + // v3: just choose the gain that gives the least "variance" in u_pid when the disturbance is simulated + // in closed loop with the given PID-model and the assumed candiate process model. + + var uPidDistance = vec.Mean(vec.Abs(vec.Diff(u_pid_adjusted))).Value; + + // v4: just choose the gain that gives the least "distance travelled" in d_est + var dEstDistance = vec.Mean(vec.Abs(vec.Diff(dEst))).Value; + + var covBtwUPidAdjustedAndDest = Math.Abs(Measures.Covariance(dEst, u_pid_adjusted, false)); + + double covarianceBtwDestAndUexternal = 0; + if (dataSet.U.GetNColumns() > 1) + { + for (int inputIdx = 0; inputIdx < dataSet.U.GetNColumns(); inputIdx++) + { + if (inputIdx == pidInputIdx) + { + continue; + } + covarianceBtwDestAndUexternal += + // Math.Abs(CorrelationCalculator.CorrelateTwoVectors(d_est, dataSet.U.GetColumn(inputIdx), dataSet.IndicesToIgnore)); + Math.Abs(Measures.Covariance(dEst, dataSet.U.GetColumn(inputIdx), false)); + // Math.Abs(CorrelationCalculator.CorrelateTwoVectors(d_est, distIdResultAlt.adjustedUnitDataSet.U.GetColumn(nonPidIdx), dataSet.IndicesToIgnore)); + } + } + unitParametersList.Add(unitParameters); - covBtwDestAndYsetList.Add(covBtwDestAndYset); + covBtwDestAndYsetList.Add(covarianceBtwDestAndYset); uPidDistanceList.Add(uPidDistance); dEstDistanceList.Add(dEstDistance); - this.covBtwDestAndUexternal.Add(covBtwDestAndUexternal); + covBtwDestAndUexternal.Add(covarianceBtwDestAndUexternal); + } /// @@ -87,16 +124,16 @@ public void Add(double gain, UnitParameters unitParameters, double covBtwDestAnd /// /// /// - /// + /// - public Tuple GetBestModel(double initalGainEstimate) + public Tuple GetBestModel() { bool doV4 = true; if (unitParametersList.Count()== 0) return new Tuple(null,""); - // calculate strenght of a minimum - strength is value between 0 and 100, higher is stronger + // calculate strength of a minimum - strength is value between 0 and 100, higher is stronger // this is a metric of how flat the objective space is, the lower the strenght, the flatter the objective function Tuple CalcStrengthOfObjectiveMinimum(double[] values) { @@ -198,8 +235,7 @@ double[] Scale(double[] v_in) return new Tuple(null, ""); } } - //modelParameters.AddWarning(UnitdentWarnings.ClosedLoopEst_GlobalSearchFailedToFindLocalMinima); - + } double[] objFun = v1; var retString = "Kp:v1(strength: " + v1_Strength.ToString("F2",CultureInfo.InvariantCulture) + ") "+ Direction(v1, min_ind_v1); ; @@ -230,5 +266,7 @@ double[] Scale(double[] v_in) return new Tuple(new UnitModel(unitParams), retString); } + + } } \ No newline at end of file diff --git a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs index 424107a..93f566e 100644 --- a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs +++ b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs @@ -108,7 +108,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters double[] u0 = SignificantDigits.Format(dataSet.U.GetRow(0), nDigits); double y0 = dataSet.Y_meas[0]; - bool isOK; + var dataSetRun1 = new UnitDataSet(dataSet); var dataSetRun2 = new UnitDataSet(dataSet); var fittingSpecs = new FittingSpecs(); @@ -140,12 +140,11 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters unitModel_step1.modelParameters.LinearGainUnc = null; idDisturbancesList.Add(distIdResult1); idUnitModelsList.Add(unitModel_step1); - isOK = ClosedLoopSim(dataSetRun1, unitModel_step1.GetModelParameters(), pidParams, distIdResult1.d_est,pidInputIdx, "run1"); if (doConsoleDebugOut) Console.WriteLine("Step1,ident: " + unitModel_step1.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture)); } - // "step2" : "global search" for linear pid-gainsgains + // "step2" : "global search" for linear pid-gains if (unitModel_step1.modelParameters.GetProcessGains() != null) { double pidProcessInputInitalGainEstimate = unitModel_step1.modelParameters.GetProcessGains()[pidInputIdx]; @@ -179,8 +178,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters // first pass(wider grid with larger grid size) var retGlobalSearch1 = GlobalSearchLinearPidGain(dataSet, pidParams, pidInputIdx, - unitModel_step1, pidProcessInputInitalGainEstimate, - min_gain, max_gain, fittingSpecs, firstPassNumIterations); + unitModel_step1, min_gain, max_gain, fittingSpecs, firstPassNumIterations); var bestUnitModel = retGlobalSearch1.Item1; if (bestUnitModel != null) { @@ -206,7 +204,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters var gainPass1 = retGlobalSearch1.Item1.modelParameters.LinearGains[pidInputIdx]; var retGlobalSearch2 = GlobalSearchLinearPidGain(dataSet, pidParams, pidInputIdx, - retGlobalSearch1.Item1, gainPass1, gainPass1 - retGlobalSearch1.Item2 * WIDTH_OF_SEARCH_PASS2, gainPass1 + retGlobalSearch1.Item2 * WIDTH_OF_SEARCH_PASS2, + retGlobalSearch1.Item1, gainPass1 - retGlobalSearch1.Item2 * WIDTH_OF_SEARCH_PASS2, gainPass1 + retGlobalSearch1.Item2 * WIDTH_OF_SEARCH_PASS2, fittingSpecs, secondPassNumIterations); bestUnitModel = retGlobalSearch2.Item1; @@ -233,12 +231,12 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters const bool doRun2 = true; const double LARGEST_TIME_CONSTANT_TO_CONSIDER_TIMEBASE_MULTIPLE = 60 + 1;//too low - // run2: do a run where it is no longer assumed that x[k-1] = y[k], + // step3 : do a run where it is no longer assumed that x[k-1] = y[k], // this run has the best chance of estimating correct time constants, but it requires a good inital guess of d - // run 2 version 1: try looking for the time constant that gives the smallest average disturbance amplitude + // version 1: try looking for the time constant that gives the smallest average disturbance amplitude - const bool doRun2Version2 = false;// TODO: implement version 2 and change over to improve Tc estimate + const int doRun2Version2 = 1;// TODO: implement version 2 and change over to improve Tc estimate if (idUnitModelsList.Count() >0) { @@ -247,7 +245,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters var unitModel = idUnitModelsList.Last(); unitModel.ID = "process"; - if (doRun2Version2) + if (doRun2Version2 == 2) { var pidModel = new PidModel(pidParams, "PID"); @@ -316,7 +314,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters // todo: determine the Tc that has the higest plant fit score Console.WriteLine(fitScores.ToString()); } - else + else // version 1 { // approach1: // "look for the time-constant that gives the disturbance that changes the least" @@ -454,7 +452,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters identUnitModel.modelParameters.Fitting.NFittingBadDataPoints = dataSet.IndicesToIgnore.Count(); // closed-loop simulation, adds U_sim and Y_sim to dataset - ClosedLoopSim(dataSet, identUnitModel.modelParameters, pidParams, disturbance,pidInputIdx); + ClosedLoopSim(dataSet, identUnitModel.modelParameters, pidParams, pidInputIdx); // round resulting parameters @@ -533,8 +531,6 @@ private static double EstimateDisturbanceLF(UnitDataSet dataSet, UnitModel unit Shared.DisablePlots(); - - return 0; } @@ -548,15 +544,13 @@ private static double EstimateDisturbanceLF(UnitDataSet dataSet, UnitModel unit /// /// /// - /// /// /// /// /// /// the best model and the step size used as a tuple - private static Tuple GlobalSearchLinearPidGain(UnitDataSet dataSet, PidParameters pidParams, int pidInputIdx, - UnitModel unitModel_run1, double pidProcessInputInitalGainEstimate, double minPidProcessGain, - double maxPidProcessGain, FittingSpecs fittingSpecs, int numberOfGlobalSearchIterations = 40) + private static Tuple GlobalSearchLinearPidGain(UnitDataSet dataSet, PidParameters pidParams, int pidInputIdx, + UnitModel unitModel_run1, double minPidProcessGain, double maxPidProcessGain, FittingSpecs fittingSpecs, int numberOfGlobalSearchIterations = 40) { var range = maxPidProcessGain - minPidProcessGain; var searchResults = new ClosedLoopGainGlobalSearchResults(); @@ -599,11 +593,9 @@ private static Tuple GlobalSearchLinearPidGain(UnitDataSet dat continue; } var dEst = curCandDisturbanceEst_SISO.d_est; - var u_pid_adjusted = curCandDisturbanceEst_SISO.adjustedUnitDataSet.U.GetColumn(pidInputIdx); - var alternativeModelParameters = curCandidateSISOModel.modelParameters.CreateCopy(); - var d_est_adjusted = dEst; + var candidateModelParameters = curCandidateSISOModel.modelParameters.CreateCopy(); // Multiple-input single-output modeling if (nGains > 1) @@ -689,7 +681,7 @@ private static Tuple GlobalSearchLinearPidGain(UnitDataSet dat } // it seems that the gains for external inputs u in model_dist are accurate if // and when pidLinProcessGain matches the "true" estimate. - d_est_adjusted = vec.Subtract(dEst, dataSet_d.Y_sim); + // d_est_adjusted = vec.Subtract(dEst, dataSet_d.Y_sim); int curModelDistInputIdx = 0; if (doIncludeYsetpointAsInput) curModelDistInputIdx = 1; @@ -711,16 +703,13 @@ private static Tuple GlobalSearchLinearPidGain(UnitDataSet dat var dataSet_altMISO = new UnitDataSet(dataSet); var distIdResultAlt_MISO = DisturbanceIdentifier.CalculateDisturbanceVector (dataSet_altMISO, alternativeMISOModel, pidInputIdx, pidParams); - var d_est_MISO = distIdResultAlt_MISO.d_est; - isOK = ClosedLoopSim(dataSet_altMISO, alternativeMISOModel.GetModelParameters(), pidParams, d_est_MISO, pidInputIdx, "run_altMISO"); - if (false) { // d_est_SISO: first pass, does not account for changes in u_externals and y_set // d_est_adj: second pass, is d_est_SISO where influence has been modelled and subtracted (this will not be perfect but should move less in repsonse to u_Ext and yset) // d_est_MISO: third pass, hopefully improves on d_est_adj- this estimate ideally does not move at all to changes in u_external and y_set - var variableList = new List { dEst, d_est_adjusted, distIdResultAlt_MISO.d_est, dataSet.Y_setpoint }; - var variableNameList = new List { "y1=d_est_SISO", "y1=d_est_adj", "y1=d_est_MISO", "y4=y_set" }; + var variableList = new List { dEst, distIdResultAlt_MISO.d_est, dataSet.Y_setpoint }; + var variableNameList = new List { "y1=d_est_SISO", "y1=d_est_MISO", "y4=y_set" }; for (int inputIdx = 0; inputIdx < nGains; inputIdx++) { if (pidInputIdx == inputIdx) @@ -733,57 +722,21 @@ private static Tuple GlobalSearchLinearPidGain(UnitDataSet dat variableList, variableNameList, dataSet_altMISO.GetTimeBase(), "ClosedLoopGlobal_d_est"); Shared.DisablePlots(); } + isOK = ClosedLoopSim(dataSet_altMISO, alternativeMISOModel.GetModelParameters(), pidParams, pidInputIdx, "run_altMISO"); + if (isOK) { - dEst = d_est_MISO; + dEst = distIdResultAlt_MISO.d_est; u_pid_adjusted = distIdResultAlt_MISO.adjustedUnitDataSet.U.GetColumn(pidInputIdx); - alternativeModelParameters = alternativeMISOModel.modelParameters.CreateCopy(); - alternativeModelParameters.Fitting = new FittingInfo(); - alternativeModelParameters.Fitting.WasAbleToIdentify = true; + candidateModelParameters = alternativeMISOModel.modelParameters.CreateCopy(); + candidateModelParameters.Fitting = new FittingInfo(); + candidateModelParameters.Fitting.WasAbleToIdentify = true; } } if (!isOK) continue; - double covarianceBtwDestAndYset = Math.Abs(CorrelationCalculator.CorrelateTwoVectors( - dEst, dataSet.Y_setpoint, dataSet.IndicesToIgnore)); - - // v3: just choose the gain that gives the least "variance" in u_pid when the disturbance is simulated - // in closed loop with the given PID-model and the assumed candiate process model. - - var uPidDistance = vec.Mean(vec.Abs(vec.Diff(u_pid_adjusted))).Value; - - // v4: just choose the gain that gives the least "distance travelled" in d_est - var dEstDistance = vec.Mean(vec.Abs(vec.Diff(dEst))).Value; // original: - // var dEstVariance = vec.Mean(vec.Diff(dEst)).Value; - //var dEstVariance = vec.Mean(vec.Abs(dEst)).Value; - //var dEstVariance = dEst.Variance(); - - // test: - var covBtwUPidAdjustedAndDest = Math.Abs(Measures.Covariance(dEst, u_pid_adjusted, false)); - - - // v5: if MISO: the disturbance should be as indpendent of the external inputs as possible - double covarianceBtwDestAndUexternal = 0; - if (dataSet.U.GetNColumns() > 1) - { - for (int inputIdx = 0; inputIdx < dataSet.U.GetNColumns(); inputIdx++) - { - if (inputIdx == pidInputIdx) - { - continue; - } - covarianceBtwDestAndUexternal += - // Math.Abs(CorrelationCalculator.CorrelateTwoVectors(d_est, dataSet.U.GetColumn(inputIdx), dataSet.IndicesToIgnore)); - Math.Abs(Measures.Covariance(dEst, dataSet.U.GetColumn(inputIdx), false)); - // Math.Abs(CorrelationCalculator.CorrelateTwoVectors(d_est, distIdResultAlt.adjustedUnitDataSet.U.GetColumn(nonPidIdx), dataSet.IndicesToIgnore)); - } - } - - - // finally,save all kpis of of this run. - searchResults.Add(curCandPidLinProcessGain, alternativeModelParameters, covarianceBtwDestAndYset, - uPidDistance, dEstDistance, covarianceBtwDestAndUexternal); + searchResults.Add(candidateModelParameters,dataSet, dEst, u_pid_adjusted, pidInputIdx); // save the time-series for debug-plotting if (doDebugPlot == true) @@ -798,8 +751,8 @@ private static Tuple GlobalSearchLinearPidGain(UnitDataSet dat uSimList.Add(u_pid_adjusted); var yProc = vec.Subtract(dataSet.Y_meas, dEst); yProcessList.Add(yProc); - - } + + } } if (doDebugPlot) @@ -822,7 +775,7 @@ double Power(double[] inSignal) double scale = Math.Max(Math.Abs(max - mean), Math.Abs(min - mean)); - return vec.Mean( vec.Abs(vec.Subtract(inSignal, vec.Mean(inSignal).Value )) ).Value/ scale; + return vec.Mean(vec.Abs(vec.Subtract(inSignal, vec.Mean(inSignal).Value))).Value / scale; } double Range(double[] inSignal) @@ -865,20 +818,20 @@ double Range(double[] inSignal) // it is acutally easier to think of disturbance d_LF and d_HF as d_y and d_u, as d_HF is really just e= ymeas-ysetp and d_u is depends on u. // this means, that d_y will not change for different guesses of the process gain or process time-constant, only d_u does that. // note that d_LF or "d_u" is basically the inverse of the y_process. - string comment = - "Kp=" + kpList.ElementAt(i).ToString("F2", CultureInfo.InvariantCulture) - + " dEstPower="+ dEstPower.ToString("F3", CultureInfo.InvariantCulture) - + " yProcPower="+yProcPower.ToString("F3", CultureInfo.InvariantCulture) - + " yMeasPower="+ yMeasPower.ToString("F3", CultureInfo.InvariantCulture) - +" uMeasPower=" + uMeasPower.ToString("F3", CultureInfo.InvariantCulture) - +" uMeanDest=" +meanDest.ToString("F3", CultureInfo.InvariantCulture); + string comment = + "Kp=" + kpList.ElementAt(i).ToString("F2", CultureInfo.InvariantCulture) + + " dEstPower=" + dEstPower.ToString("F3", CultureInfo.InvariantCulture) + + " yProcPower=" + yProcPower.ToString("F3", CultureInfo.InvariantCulture) + + " yMeasPower=" + yMeasPower.ToString("F3", CultureInfo.InvariantCulture) + + " uMeasPower=" + uMeasPower.ToString("F3", CultureInfo.InvariantCulture) + + " uMeanDest=" + meanDest.ToString("F3", CultureInfo.InvariantCulture); Shared.EnablePlots(); - Plot.FromList(new List { dEstList.ElementAt(i), vec.Subtract(yProcessList.ElementAt(i), vec.Mean(yProcessList.ElementAt(i)).Value), + Plot.FromList(new List { dEstList.ElementAt(i), vec.Subtract(yProcessList.ElementAt(i), vec.Mean(yProcessList.ElementAt(i)).Value), vec.Subtract(dataSet.Y_meas, vec.Mean(dataSet.Y_meas).Value), vec.Subtract(dataSet.U.GetColumn(0),vec.Mean(dataSet.U.GetColumn(0)).Value), dataSet.U.GetColumn(0) }, - new List { "y1=d_est", "y1=y_process", "y1=y_meas","y1=u_pidDetrend","y3=u_pid" }, dataSet.Times,comment); + new List { "y1=d_est", "y1=y_process", "y1=y_meas", "y1=u_pidDetrend", "y3=u_pid" }, dataSet.Times, comment); Shared.DisablePlots(); Thread.Sleep(100); } @@ -886,7 +839,7 @@ double Range(double[] inSignal) // finally, select the best model from all the tested models by looking across all the kpis stored - (UnitModel bestUnitModel,string gsMiniumType) = searchResults.GetBestModel(pidProcessInputInitalGainEstimate); + (var bestUnitModel, var gsMiniumType) = searchResults.GetBestModel(); if (bestUnitModel != null) { if (bestUnitModel.modelParameters != null) @@ -905,9 +858,11 @@ double Range(double[] inSignal) bestUnitModel.modelParameters.Fitting.WasAbleToIdentify = true; bestUnitModel.modelParameters.Fitting.SolverID = gsMiniumType; } - return new Tuple(bestUnitModel, gainStepSize); + return new Tuple(bestUnitModel, gainStepSize); } + + /// /// Guess the sign of the pid-input /// @@ -1082,7 +1037,7 @@ private static Tuple EstimateSISOdisturbanceForP var d_est = candidateModelDisturbance.d_est; var isOK = ClosedLoopSim - (dataSet_SISO, candidateModel.GetModelParameters(), pidParams, d_est, pidInputIdx_NewSystem, "run_altSISO"); + (dataSet_SISO, candidateModel.GetModelParameters(), pidParams, pidInputIdx_NewSystem, "run_altSISO"); if (isOK) return new Tuple(candidateModel, candidateModelDisturbance); @@ -1100,19 +1055,17 @@ private static Tuple EstimateSISOdisturbanceForP /// unit dataset for simulation /// paramters or UnitModel /// parameters of PidModel - /// disturbance vector /// index of the pid input in the process model /// optional name used for plotting /// public static bool ClosedLoopSim(UnitDataSet unitData, UnitParameters modelParams, PidParameters pidParams, - double[] disturbance,int pidInputIdx, string name="") + int pidInputIdx, string name="") { if (pidParams == null) { return false; } - // unitData.D = disturbance; var pidModel = new PidModel(pidParams,"PidModel"); var unitModel = new UnitModel(modelParams, "ProcModel"); (var plantSim, var inputDataSet) = PlantSimulator.CreateFeedbackLoopWithEstimatedDisturbance(unitData, pidModel, diff --git a/Dynamic/Identification/DisturbanceIdentifier.cs b/Dynamic/Identification/DisturbanceIdentifier.cs index c56627d..7ad0f14 100644 --- a/Dynamic/Identification/DisturbanceIdentifier.cs +++ b/Dynamic/Identification/DisturbanceIdentifier.cs @@ -332,9 +332,6 @@ public static DisturbanceIdResult CalculateDisturbanceVector(UnitDataSet unitDat // Plot.FromList(new List { d_est, d_est2, d_HF }, new List { "y1=d1", "y1=d2","y1=d_HF" }, unitDataSet.Times, "testplot"); // Shared.DisablePlots(); - - - result.d_est = d_est; result.d_LF = d_LF; result.d_HF = d_HF; diff --git a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_DynamicSISO.cs b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_DynamicSISO.cs index a5aa881..342d57f 100644 --- a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_DynamicSISO.cs +++ b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_DynamicSISO.cs @@ -117,7 +117,7 @@ public void SinusDisturbance(double distSinusAmplitude) public void RandomWalkDisturbance(double procGain, double distAmplitude, double gainPrecisionPrc) { - int seed = 150; + int seed = 10; int N = 2000; @@ -177,7 +177,7 @@ public void LongStepDist_EstimatesOk(double stepAmplitude,double procGainAllowed // for some reason, this test also does not work with "run all", but works fine when run individually? - [Test,Explicit,NonParallelizable] + [TestCase( Category = "NotWorking_AcceptanceTest"),Explicit] public void StepAtStartOfDataset_IsExcludedFromAnalysis() { double stepAmplitude = 1; diff --git a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs index 1c69b64..188305e 100644 --- a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs +++ b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs @@ -1,5 +1,6 @@ using System; using System.Collections.Generic; +using System.ComponentModel; using System.Linq; using System.Text; @@ -187,7 +188,7 @@ public void FlatData_DoesNotCrash() trueDisturbance,false,false); } - [Test] + [TestCase(Category = "NotWorking_AcceptanceTest")] public void BadDataPoints_IsExcludedFromAnalysis() { double stepAmplitude = 1; diff --git a/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs b/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs index c6b5b03..a168ab9 100644 --- a/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs +++ b/TimeSeriesAnalysis.Tests/Tests/DisturbanceEstimatorTests.cs @@ -59,13 +59,14 @@ public void CommonPlotAndAsserts(UnitDataSet pidDataSet, double[] d_est, double[ Assert.IsTrue(d_est != null); string caseId = TestContext.CurrentContext.Test.Name.Replace("(", "_"). Replace(")", "_").Replace(",", "_") + "y"; - - Plot.FromList(new List{ pidDataSet.Y_meas, pidDataSet.Y_setpoint, + if (false) + { + Plot.FromList(new List{ pidDataSet.Y_meas, pidDataSet.Y_setpoint, pidDataSet.U.GetColumn(0), d_est, trueDisturbance }, - new List { "y1=y meas", "y1=y set", "y2=u(right)", "y3=est disturbance", "y3=true disturbance" }, - pidDataSet.GetTimeBase(), caseId); - + new List { "y1=y meas", "y1=y set", "y2=u(right)", "y3=est disturbance", "y3=true disturbance" }, + pidDataSet.GetTimeBase(), caseId); + } Assert.IsTrue(vec.Mean(vec.Abs(vec.Subtract(trueDisturbance, d_est))) < distTrueAmplitude / 10,"true disturbance and actual disturbance too far apart"); } @@ -181,14 +182,6 @@ public void DisturbanceTestUsingPlantSimulator(UnitModel processModel, double[] { var pidModel1 = new PidModel(pidParameters1, "PID1"); - // important to add a disturbance signal with naming convention to the output of the process, to signal that the disturbance is to be estimated! - - // begin testing - // var dynamicModelParameters_new = dynamicModelParameters.CreateCopy(); - // dynamicModelParameters_new.TimeConstant_s = 20; - //end testing - - //var processModel3 = new UnitModel(dynamicModelParameters_new, "Proc1"); var processModel3 = new UnitModel(dynamicModelParameters, "Proc1"); var distSignal = SignalNamer.EstDisturbance(processModel3); diff --git a/docs/articles/sysid_disturbance.md b/docs/articles/sysid_disturbance.md index 99e5d55..26b1870 100644 --- a/docs/articles/sysid_disturbance.md +++ b/docs/articles/sysid_disturbance.md @@ -208,7 +208,7 @@ Because no process dynamics are assumed yet, ``d_est(t)`` at this stage will inc has dynamics. ``d_est(t)`` will be attempted improved in subsquent steps. -### Guessing the sign of the process gain +#### Guessing the sign of the process gain Methods for open-loop estimation when applied naively to closed loop time-series often estimate process gain with incorrect sign. @@ -283,19 +283,10 @@ expressed as: Some obervations: - this method is heurisitic -- the objective has a minumum +- the objective usually has a minimum, but not always (such as if the disturbance is a perfect sinus in unit tests) - the objective space is fairly flat, the minum has a fairly low ``strength'', i.e neighborhing process gains have almost equally low objectives - the objective space seems to be more concave ("stronger" i.e. more significant minimums) when the process gain is higher. -It could be an idea to count the number of zero crossing of e in the dataset, if it is very low, that may be an indication that the process gain is low. - -*It could be important to look at the internal output of the process model.* - -*It could be that it is possible to have a very good estiamte of the distrubance d, yet the process gain can still be of by 40 percent. -This has been observed in some random walk unit tests where the process gain is small. This could be because the d_est in this cae is mostly -determined by ``d_HF```? This could be the clue to determining the uncertainty of the process gain, looking at how much the disturbance vector -estimate is sensitive to the process gain values?* - There are essentially two ways of calculating the disturbance 1. By subtracting the modelled ``y_{proc}(\hat{u})`` from ``\bar{y}'' : ``d_est = \bar{y} - y_{proc}(\hat{u})`` @@ -312,7 +303,6 @@ Combining the two above means that The smaller the process gain is the more d_est is similar to d_LF. - It is possible to plot the solution space of the d_est for different Kp, and in periods where there is small changes in the integral term of the pid-controller, the disturbance looks quite similar for different Kp. So while in some periods vary with a factor 10 when Kp varies with a factor 10, in other peridos it just varies with a factor 2. @@ -321,8 +311,6 @@ when Kp varies with a factor 10, in other peridos it just varies with a factor 2 - - ### Estimating upper-and lower bounds on the process gain It may be that the algorithm has a better accuracy when there are setpoint changes in the dataset. @@ -336,7 +324,10 @@ Knowing upper- and lower bounds could also be a useful guide for subsequent manu -### Determing the dynamic parameters of the process + + + +### Step3: Determing the dynamic parameters of the process If the process is actually dynamic yet is modeled as static, then the above methodology will result into un-modeled transients bleeding into the estimated disturbance, where they will appear as