From 6c715c3cbe0adada45552e85d2caaef8833e9f3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Steinar=20Elgs=C3=A6ter?= Date: Fri, 13 Dec 2024 15:17:13 +0100 Subject: [PATCH] refactoring CLUI to allow multiple passes of algorithm --- .../ClosedLoopUnitIdentifier.cs | 366 +++++++++--------- .../Tests/ClosedLoopIdentifierTester_MISO.cs | 22 +- .../ClosedLoopIdentifierTester_StaticSISO.cs | 56 +-- 3 files changed, 214 insertions(+), 230 deletions(-) diff --git a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs index 93f566e..d60bf66 100644 --- a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs +++ b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs @@ -110,90 +110,85 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters double y0 = dataSet.Y_meas[0]; var dataSetRun1 = new UnitDataSet(dataSet); - var dataSetRun2 = new UnitDataSet(dataSet); + var dataSetStep3 = new UnitDataSet(dataSet); var fittingSpecs = new FittingSpecs(); fittingSpecs.u0 = u0; var GSdescription = ""; - { - // ---------------- - // run1: no process model assumed, let disturbance estimator guesstimate a pid-process gain, - // to give afirst estimate of the disturbance - var unitModel_step1 = ModelFreeEstimateClosedLoopProcessGain(dataSetRun1, pidParams, pidInputIdx); - var distIdResult1 = DisturbanceIdentifier.CalculateDisturbanceVector(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); + // ---------------- + // step 1: no process model assumed, let disturbance estimator guesstimate a pid-process gain, + // to give afirst estimate of the disturbance + var unitModel_step1 = ModelFreeEstimateClosedLoopProcessGain(dataSetRun1, pidParams, pidInputIdx); + var distIdResult1 = DisturbanceIdentifier.CalculateDisturbanceVector(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); + if (doConsoleDebugOut) + ConsoleDebugOut(unitModel_step1, "Step 1"); + idUnitModelsList.Add(unitModel_step1); + + // step1, MISO: ident (attempt to identify any other inputs) + if (isMISO) + { + dataSetRun1.D = distIdResult1.d_est; + unitModel_step1 = UnitIdentifier.IdentifyLinearAndStatic(ref dataSetRun1, fittingSpecs, doTimeDelayEstOnRun1); + unitModel_step1.modelParameters.LinearGainUnc = null; + idDisturbancesList.Add(distIdResult1); + idUnitModelsList.Add(unitModel_step1); if (doConsoleDebugOut) - Console.WriteLine("Step1: " + unitModel_step1.GetModelParameters().LinearGains.ElementAt(pidInputIdx).ToString("F3", CultureInfo.InvariantCulture)); + ConsoleDebugOut(unitModel_step1, "Step 1,MISO"); + } - idUnitModelsList.Add(unitModel_step1); + // var KPest = EstimateDisturbanceLF(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); - // var KPest = EstimateDisturbanceLF(dataSetRun1, unitModel_step1, pidInputIdx, pidParams); - // if (doConsoleDebugOut) - // Console.WriteLine("experimental: " + KPest); - - // run1: ident (attempt to identify any other inputs) - if (isMISO) - { - dataSetRun1.D = distIdResult1.d_est; - unitModel_step1 = UnitIdentifier.IdentifyLinearAndStatic(ref dataSetRun1, fittingSpecs, doTimeDelayEstOnRun1); - unitModel_step1.modelParameters.LinearGainUnc = null; - idDisturbancesList.Add(distIdResult1); - idUnitModelsList.Add(unitModel_step1); - if (doConsoleDebugOut) - Console.WriteLine("Step1,ident: " + unitModel_step1.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture)); - } + int MAX_RUNS = 1; + for (int runIdx = 1;runIdx<= MAX_RUNS; runIdx++) + { // "step2" : "global search" for linear pid-gains - if (unitModel_step1.modelParameters.GetProcessGains() != null) + if (idUnitModelsList.Last().modelParameters.GetProcessGains() != null) { - double pidProcessInputInitalGainEstimate = unitModel_step1.modelParameters.GetProcessGains()[pidInputIdx]; + double pidProcessInputInitalGainEstimate = idUnitModelsList.Last().modelParameters.GetProcessGains()[pidInputIdx]; var gainAndCorrDict = new Dictionary(); // // looking to find the process gain that "decouples" d_est from Y_setpoint as much as possible. - // // or antoher way to look at ti is that the output U with setpoint effects removed should be as decoupled // form Y_setpoint. - double max_gain, min_gain; - // todo: considre improving these bounds? + // todo: consider improving these bounds? if (pidProcessInputInitalGainEstimate > 0) { max_gain = pidProcessInputInitalGainEstimate * initalGuessFactor_higherbound; min_gain = pidProcessInputInitalGainEstimate * 1 / initalGuessFactor_higherbound; // mostly works, but not for sinus disturbances - //min_gain = 0; } else { min_gain = pidProcessInputInitalGainEstimate * initalGuessFactor_higherbound; max_gain = pidProcessInputInitalGainEstimate * 1 / initalGuessFactor_higherbound; - // max_gain = 0; } if (doConsoleDebugOut) { - Console.WriteLine("Step2,GS : " + min_gain.ToString("F3", CultureInfo.InvariantCulture) + " to " + max_gain.ToString("F3", CultureInfo.InvariantCulture)); + Console.WriteLine("run"+ runIdx + " Step2 " + "bounds: " + min_gain.ToString("F3", CultureInfo.InvariantCulture) + " to " + max_gain.ToString("F3", CultureInfo.InvariantCulture)); } // first pass(wider grid with larger grid size) var retGlobalSearch1 = GlobalSearchLinearPidGain(dataSet, pidParams, pidInputIdx, - unitModel_step1, min_gain, max_gain, fittingSpecs, firstPassNumIterations); + idUnitModelsList.Last(), min_gain, max_gain, fittingSpecs, firstPassNumIterations); var bestUnitModel = retGlobalSearch1.Item1; if (bestUnitModel != null) { GSdescription = bestUnitModel.GetModelParameters().Fitting.SolverID; wasGainGlobalSearchDone = true; } - + if (doConsoleDebugOut && retGlobalSearch1.Item1 != null) { - - Console.WriteLine("Step2,GS1: " + retGlobalSearch1.Item1.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture)); + ConsoleDebugOut( bestUnitModel, "run" + runIdx + "Step 2"); } else - Console.WriteLine("Step2,GS1: FAILED"); + Console.WriteLine("run" + runIdx + "Step2,GS1: FAILED"); if (bestUnitModel != null) { @@ -201,15 +196,14 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters if (retGlobalSearch1.Item1.modelParameters.Fitting.WasAbleToIdentify && secondPassNumIterations > 0) { const int WIDTH_OF_SEARCH_PASS2 = 3; - var gainPass1 = retGlobalSearch1.Item1.modelParameters.LinearGains[pidInputIdx]; var retGlobalSearch2 = GlobalSearchLinearPidGain(dataSet, pidParams, pidInputIdx, retGlobalSearch1.Item1, gainPass1 - retGlobalSearch1.Item2 * WIDTH_OF_SEARCH_PASS2, gainPass1 + retGlobalSearch1.Item2 * WIDTH_OF_SEARCH_PASS2, fittingSpecs, secondPassNumIterations); bestUnitModel = retGlobalSearch2.Item1; - + if (doConsoleDebugOut) - Console.WriteLine("Step2,GS2: " + retGlobalSearch2.Item1.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture)); + ConsoleDebugOut(bestUnitModel, "run" + runIdx + "Step 2,GS2"); } } // add the "best" model to be used in the next model run @@ -218,159 +212,162 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters idUnitModelsList.Add(bestUnitModel); } } - } + - // ---------------- - // - 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 - // time delays are way off if the process is excited mainly by disturbances. - // - the reason that we cannot do run2 immediately, is that that formulation - // does not appear to give a solution if the guess disturbance vector is bad. + // ---------------- + // - 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 + // time delays are way off if the process is excited mainly by disturbances. + // - the reason that we cannot do run2 immediately, is that that formulation + // does not appear to give a solution if the guess disturbance vector is bad. - const bool doRun2 = true; - const double LARGEST_TIME_CONSTANT_TO_CONSIDER_TIMEBASE_MULTIPLE = 60 + 1;//too low + const bool doStep3 = true; + const double LARGEST_TIME_CONSTANT_TO_CONSIDER_TIMEBASE_MULTIPLE = 60 + 1;//too low - // 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 + // 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 - // 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 int doRun2Version2 = 1;// TODO: implement version 2 and change over to improve Tc estimate + const int doStep3Version2 = 1;// TODO: implement version 2 and change over to improve Tc estimate - if (idUnitModelsList.Count() >0) - { - if (doRun2 && idUnitModelsList.Last() != null) + if (idUnitModelsList.Count() > 0) { - var unitModel = idUnitModelsList.Last(); - unitModel.ID = "process"; - - if (doRun2Version2 == 2) + if (doStep3 && idUnitModelsList.Last() != null) { - var pidModel = new PidModel(pidParams, "PID"); - - var fitScores = new List(); - var Tcs = new List(); - - var TcListTest = new List() { 0, 5, 10, 15 };//for debugging, generalize - var modelParams = unitModel.GetModelParameters(); - - // TODO: need to estimate disturbance and add it to the simulation in each case. - var runCounter = 0; - var d_estList = new List(); - var d_estDesc = new List(); - var y_simList = new List(); - var y_simDesc = new List(); - var u_simList = new List(); - var u_simDesc = new List(); - DateTime[] dateTimes = null; - bool doDebugPlot = true; - - (var plantSim, var inputData) = PlantSimulator.CreateFeedbackLoopWithEstimatedDisturbance(dataSetRun2, pidModel, unitModel, pidInputIdx); - if (doDebugPlot) - { - y_simList.Add(inputData.GetValues(unitModel.ID, SignalType.Output_Y)); - y_simDesc.Add("y1= y_meas"); - u_simList.Add(inputData.GetValues(pidModel.ID, SignalType.PID_U)); - u_simDesc.Add("y3= u_meas"); - } + var unitModel = idUnitModelsList.Last(); + unitModel.ID = "process"; - foreach (var Tc in TcListTest) + if (doStep3Version2 == 2) { - modelParams.TimeConstant_s = Tc; - ((UnitModel)plantSim.modelDict[unitModel.ID]).SetModelParameters(modelParams); + var pidModel = new PidModel(pidParams, "PID"); + + var fitScores = new List(); + var Tcs = new List(); + + var TcListTest = new List() { 0, 5, 10, 15 };//for debugging, generalize + var modelParams = unitModel.GetModelParameters(); + + // TODO: need to estimate disturbance and add it to the simulation in each case. + var runCounter = 0; + var d_estList = new List(); + var d_estDesc = new List(); + var y_simList = new List(); + var y_simDesc = new List(); + var u_simList = new List(); + var u_simDesc = new List(); + DateTime[] dateTimes = null; + bool doDebugPlot = true; + + (var plantSim, var inputData) = PlantSimulator.CreateFeedbackLoopWithEstimatedDisturbance(dataSetStep3, pidModel, unitModel, pidInputIdx); + if (doDebugPlot) + { + y_simList.Add(inputData.GetValues(unitModel.ID, SignalType.Output_Y)); + y_simDesc.Add("y1= y_meas"); + u_simList.Add(inputData.GetValues(pidModel.ID, SignalType.PID_U)); + u_simDesc.Add("y3= u_meas"); + } - var inputDataCoSim = new TimeSeriesDataSet(inputData); - // 2. co-simulate with disturbance from 1. - bool isCoSimOK = plantSim.Simulate(inputDataCoSim, out var simData); - if (isCoSimOK) + foreach (var Tc in TcListTest) { - fitScores.Add(plantSim.PlantFitScore); // issue: these are now always 100%??? - Tcs.Add(Tc); + modelParams.TimeConstant_s = Tc; + ((UnitModel)plantSim.modelDict[unitModel.ID]).SetModelParameters(modelParams); + + var inputDataCoSim = new TimeSeriesDataSet(inputData); + // 2. co-simulate with disturbance from 1. + bool isCoSimOK = plantSim.Simulate(inputDataCoSim, out var simData); + if (isCoSimOK) + { + fitScores.Add(plantSim.PlantFitScore); // issue: these are now always 100%??? + Tcs.Add(Tc); + } + // debugging plot + if (doDebugPlot) + { + y_simList.Add(simData.GetValues(unitModel.ID, SignalType.Output_Y)); + y_simDesc.Add("y1= y_sim" + runCounter); + u_simList.Add(simData.GetValues(pidModel.ID, SignalType.PID_U)); + u_simDesc.Add("y3= u_sim" + runCounter); + d_estList.Add(simData.GetValues(unitModel.ID, SignalType.Disturbance_D)); + d_estDesc.Add("y2= d_est" + runCounter); + dateTimes = inputDataCoSim.GetTimeStamps(); + } + runCounter++; } - // debugging plot if (doDebugPlot) { - y_simList.Add(simData.GetValues(unitModel.ID, SignalType.Output_Y)); - y_simDesc.Add("y1= y_sim" + runCounter); - u_simList.Add(simData.GetValues(pidModel.ID, SignalType.PID_U)); - u_simDesc.Add("y3= u_sim" + runCounter); - d_estList.Add(simData.GetValues(unitModel.ID, SignalType.Disturbance_D)); - d_estDesc.Add("y2= d_est" + runCounter); - dateTimes = inputDataCoSim.GetTimeStamps(); + y_simList.AddRange(u_simList); + y_simList.AddRange(d_estList); + y_simDesc.AddRange(u_simDesc); + y_simDesc.AddRange(d_estDesc); + Shared.EnablePlots(); + Plot.FromList(y_simList, y_simDesc, dateTimes, "debug"); + Shared.DisablePlots(); } - runCounter++; + // todo: determine the Tc that has the higest plant fit score + Console.WriteLine(fitScores.ToString()); } - if (doDebugPlot) + else // version 1 { - y_simList.AddRange(u_simList); - y_simList.AddRange(d_estList); - y_simDesc.AddRange(u_simDesc); - y_simDesc.AddRange(d_estDesc); - Shared.EnablePlots(); - Plot.FromList(y_simList, y_simDesc, dateTimes, "debug"); - Shared.DisablePlots(); - } - // todo: determine the Tc that has the higest plant fit score - Console.WriteLine(fitScores.ToString()); - } - else // version 1 - { - // approach1: - // "look for the time-constant that gives the disturbance that changes the least" - // - var distIdResult2 = DisturbanceIdentifier.CalculateDisturbanceVector - (dataSetRun2, unitModel, pidInputIdx, pidParams); - var estDisturbances = new List(); - var distDevs = new List(); - - estDisturbances.Add(distIdResult2.d_est); - var timeBase = dataSetRun2.GetTimeBase(); - double candiateTc_s = 0; - bool curDevIsDecreasing = true; - double firstDev = vec.Sum(vec.Abs(vec.Diff(distIdResult2.d_est))).Value; - distDevs.Add(firstDev); - while (candiateTc_s < LARGEST_TIME_CONSTANT_TO_CONSIDER_TIMEBASE_MULTIPLE * timeBase && curDevIsDecreasing) - { - candiateTc_s += timeBase; - var newParams = unitModel.GetModelParameters().CreateCopy(); - newParams.TimeConstant_s = candiateTc_s; - var newModel = new UnitModel(newParams); - var distIdResult_Test = DisturbanceIdentifier.CalculateDisturbanceVector - (dataSetRun2, newModel, pidInputIdx, pidParams); - estDisturbances.Add(distIdResult_Test.d_est); - var curDev = vec.Sum(vec.Abs(vec.Diff(distIdResult_Test.d_est))).Value; - if (curDev < distDevs.Last()) - curDevIsDecreasing = true; - else - curDevIsDecreasing = false; - distDevs.Add(curDev); - } + // approach1: + // "look for the time-constant that gives the disturbance that changes the least" + // + var distIdResult2 = DisturbanceIdentifier.CalculateDisturbanceVector + (dataSetStep3, unitModel, pidInputIdx, pidParams); + var estDisturbances = new List(); + var distDevs = new List(); + + estDisturbances.Add(distIdResult2.d_est); + var timeBase = dataSetStep3.GetTimeBase(); + double candiateTc_s = 0; + bool curDevIsDecreasing = true; + double firstDev = vec.Sum(vec.Abs(vec.Diff(distIdResult2.d_est))).Value; + distDevs.Add(firstDev); + while (candiateTc_s < LARGEST_TIME_CONSTANT_TO_CONSIDER_TIMEBASE_MULTIPLE * timeBase && curDevIsDecreasing) + { + candiateTc_s += timeBase; + var newParams = unitModel.GetModelParameters().CreateCopy(); + newParams.TimeConstant_s = candiateTc_s; + var newModel = new UnitModel(newParams); + var distIdResult_Test = DisturbanceIdentifier.CalculateDisturbanceVector + (dataSetStep3, newModel, pidInputIdx, pidParams); + estDisturbances.Add(distIdResult_Test.d_est); + + var curDev = vec.Mean(vec.Abs(vec.Diff(distIdResult_Test.d_est))).Value; + + if (curDev < distDevs.Last()) + curDevIsDecreasing = true; + else + curDevIsDecreasing = false; + distDevs.Add(curDev); + } - if (candiateTc_s == LARGEST_TIME_CONSTANT_TO_CONSIDER_TIMEBASE_MULTIPLE) - { - // you may get here when the disturbance is continous and noisy - Console.WriteLine("warning: ClosedLoopIdentifierFailedToFindAUniqueProcessTc"); - } + if (candiateTc_s == LARGEST_TIME_CONSTANT_TO_CONSIDER_TIMEBASE_MULTIPLE) + { + // you may get here when the disturbance is continous and noisy + Console.WriteLine("warning: ClosedLoopIdentifierFailedToFindAUniqueProcessTc"); + } + + // TODO: it would be possible to divide the time-constant into a time-delay and a time constant + if (candiateTc_s > 0) + { + candiateTc_s -= timeBase; + } + var step3params = unitModel.GetModelParameters().CreateCopy(); + step3params.TimeConstant_s = candiateTc_s; + var step3Model = new UnitModel(step3params); + idUnitModelsList.Add(step3Model); + var distIdResult_step4 = DisturbanceIdentifier.CalculateDisturbanceVector + (dataSetStep3, step3Model, pidInputIdx, pidParams); + idDisturbancesList.Add(distIdResult_step4); + if (doConsoleDebugOut) + ConsoleDebugOut(step3Model, "run" + runIdx + "Step 3"); - // TODO: it would be possible to divide the time-constant into a time-delay and a time constant - if (candiateTc_s > 0) - { - candiateTc_s -= timeBase; } - var step2params = unitModel.GetModelParameters().CreateCopy(); - step2params.TimeConstant_s = candiateTc_s; - var step2Model = new UnitModel(step2params); - idUnitModelsList.Add(step2Model); - var distIdResult_step4 = DisturbanceIdentifier.CalculateDisturbanceVector - (dataSetRun2, step2Model, pidInputIdx, pidParams); - idDisturbancesList.Add(distIdResult_step4); - if (doConsoleDebugOut) - Console.WriteLine("Run2 " + step2Model.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture)); } } } - // debugging plots, should normally be off @@ -382,7 +379,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters Console.WriteLine("run2"); Console.WriteLine(idUnitModelsList[1].ToString()); Plot.FromList( - new List { + new List { idDisturbancesList[0].d_est, idDisturbancesList[1].d_est, @@ -392,14 +389,14 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters idDisturbancesList[0].d_LF, idDisturbancesList[1].d_LF, }, - new List {"y1=d_run1", "y1=d_run2", + new List {"y1=d_run1", "y1=d_run2", "y3=dHF_run1", "y3=dHF_run2", "y3=dLF_run1", "y3=dLF_run2", }, dataSet.GetTimeBase(), "doDebuggingPlot_disturbanceHLandLF"); dataSet.D = null; - (var isOk1,var sim1results) = PlantSimulator.SimulateSingle(dataSet, idUnitModelsList[0], false); + (var isOk1, var sim1results) = PlantSimulator.SimulateSingle(dataSet, idUnitModelsList[0], false); (var isOk2, var sim2results) = PlantSimulator.SimulateSingle(dataSet, idUnitModelsList[1], false); (var isOk3, var sim3results) = PlantSimulator.SimulateSingle(dataSet, idUnitModelsList[2], false); @@ -414,7 +411,7 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters new List {"y1=y_meas","y1=xd_run1", "y1=xd_run2", "y3=x_run1","y3=x_run2"}, dataSet.GetTimeBase(), "doDebuggingPlot_states"); - Shared.DisablePlots(); + Shared.DisablePlots(); } if (idDisturbancesList.Count > 0) @@ -444,13 +441,13 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters if (wasGainGlobalSearchDone) { - identUnitModel.modelParameters.Fitting.SolverID = "ClosedLoop/w gain global search "+ GSdescription; + identUnitModel.modelParameters.Fitting.SolverID = "ClosedLoop/w gain global search " + GSdescription; } else identUnitModel.modelParameters.Fitting.SolverID = "ClosedLoop heuritic (global search no minimum found)"; identUnitModel.modelParameters.Fitting.NFittingTotalDataPoints = dataSet.GetNumDataPoints(); identUnitModel.modelParameters.Fitting.NFittingBadDataPoints = dataSet.IndicesToIgnore.Count(); - + // closed-loop simulation, adds U_sim and Y_sim to dataset ClosedLoopSim(dataSet, identUnitModel.modelParameters, pidParams, pidInputIdx); @@ -479,10 +476,10 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters return (retModel, null); } - + } - private static double EstimateDisturbanceLF(UnitDataSet dataSet, UnitModel unitModel, int pidInputIdx, PidParameters pidParams) + private static double EstimateDisturbanceLF(UnitDataSet dataSet, UnitModel unitModel, int pidInputIdx, PidParameters pidParams) { var vec = new Vec(dataSet.BadDataID); @@ -509,14 +506,14 @@ private static double EstimateDisturbanceLF(UnitDataSet dataSet, UnitModel unit umInternal.SetModelParameters(unitParams); (var isOk, var y_proc) = PlantSimulator.SimulateSingle(dataSet, umInternal, false); var d_LF = vec.Multiply(vec.Subtract(y_proc, y_proc[0]), -1); - var d_est1 = vec.Add(d_HF, d_LF); + var d_est1 = vec.Add(d_HF, d_LF); var d_est2 = vec.Subtract(dataSet.Y_meas, y_proc); dList.Add(d_est2); - dNameList.Add("y1=destKp"+Kp.ToString("F2",CultureInfo.InvariantCulture)); - // var dLF_energy = vec.Mean(d_LF).Value; + dNameList.Add("y1=destKp" + Kp.ToString("F2", CultureInfo.InvariantCulture)); + // var dLF_energy = vec.Mean(d_LF).Value; //Debug.WriteLine("Kp="+Kp.ToString("F2", CultureInfo.InvariantCulture) + "dLFenergy:"+ dLF_energy.ToString("F2", CultureInfo.InvariantCulture) - // + "dHFenergy" +dHF_energy.ToString("F2", CultureInfo.InvariantCulture)) ; + // + "dHFenergy" +dHF_energy.ToString("F2", CultureInfo.InvariantCulture)) ; v2List.Add(y_proc); v2NameList.Add("y3=yProcKp" + Kp.ToString("F2", CultureInfo.InvariantCulture)); @@ -529,11 +526,16 @@ private static double EstimateDisturbanceLF(UnitDataSet dataSet, UnitModel unit Shared.EnablePlots(); Plot.FromList(dList, dNameList, dataSet.GetTimeBase(), "exp_Kp"); Shared.DisablePlots(); - - return 0; } + private static void ConsoleDebugOut(UnitModel unitModel, string description) + { + Console.WriteLine(description + " Kp:" + unitModel.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture) + + " Tc:" + unitModel.GetModelParameters().TimeConstant_s.ToString("F1", CultureInfo.InvariantCulture)); + } + + /// /// tries N=numberOfGlobalSearchIterations process gains of the unitModel in the range[minPidProcessGain,maxPidProcessGain], using a closed-loop simulation /// that includes the PID-model with paramters pidParams that acts on input pidInputIdx of the unit model. diff --git a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_MISO.cs b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_MISO.cs index 756fe3e..4844c80 100644 --- a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_MISO.cs +++ b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_MISO.cs @@ -16,6 +16,7 @@ namespace TimeSeriesAnalysis.Test.DisturbanceID class ClosedLoopIdentifierTester_MISO { const bool doPlot = true; + const bool isStatic = true; UnitParameters staticModelParameters = new UnitParameters { @@ -52,7 +53,7 @@ public void SetUp() } public void CommonPlotAndAsserts(UnitDataSet pidDataSet, double[] estDisturbance, double[] trueDisturbance, - UnitModel identifiedModel, UnitModel trueModel, double maxAllowedGainOffsetPrc, double maxAllowedMeanDisturbanceOffsetPrc = 30) + UnitModel identifiedModel, UnitModel trueModel, double maxAllowedGainOffsetPrc, double maxAllowedMeanDisturbanceOffsetPrc = 30, bool isStatic = true) { Vec vec = new Vec(); @@ -89,10 +90,13 @@ public void CommonPlotAndAsserts(UnitDataSet pidDataSet, double[] estDisturbance } // time constant: - double estTimeConstant_s = identifiedModel.modelParameters.TimeConstant_s; - double trueTimeConstant_s = trueModel.modelParameters.TimeConstant_s; - double timeConstant_tol_s = 1; - Assert.IsTrue(Math.Abs(estTimeConstant_s - trueTimeConstant_s) < timeConstant_tol_s, "est time constant " + estTimeConstant_s + " too far off " + trueTimeConstant_s); + if (!isStatic) + { + double estTimeConstant_s = identifiedModel.modelParameters.TimeConstant_s; + double trueTimeConstant_s = trueModel.modelParameters.TimeConstant_s; + double timeConstant_tol_s = 1; + Assert.IsTrue(Math.Abs(estTimeConstant_s - trueTimeConstant_s) < timeConstant_tol_s, "est time constant " + estTimeConstant_s + " too far off " + trueTimeConstant_s); + } } [TestCase(0,false)] @@ -113,7 +117,7 @@ public void StaticMISO_SetpointChanges_no_disturbance_detectsProcessOk(int pidIn trueParamters.LinearGains = new double[] { oldGains[1], oldGains[0], oldGains[2] }; } GenericMISODisturbanceTest(new UnitModel(trueParamters, "StaticProcess"), trueDisturbance, externalU1,externalU2, - doNegative, true,yset, pidInputIdx); + doNegative, true,yset, pidInputIdx,10,false, isStatic); } [TestCase(0, false, false), NonParallelizable] @@ -159,7 +163,7 @@ public void StaticMISO_SetpointChanges_WITH_disturbance_detectsProcessOk(int pid var externalU1 = TimeSeriesCreator.Step(N / 8, N, 15, 20); var yset = TimeSeriesCreator.Step(N * 3 / 8, N, 20, 18); GenericMISODisturbanceTest(new UnitModel(trueParamters, "StaticProcess"), trueDisturbance, externalU1, externalU2, - doNegative, true, yset, pidInputIdx); + doNegative, true, yset, pidInputIdx,10,false,isStatic); } // be aware that adding any sort of dynamics to the "true" model here seems to destroy the // model estimate. @@ -240,7 +244,7 @@ public void DynamicMISO_externalUchanges_NoDisturbance_NOsetpointchange_DetectsP public void GenericMISODisturbanceTest (UnitModel trueProcessModel, double[] trueDisturbance, double[] externalU1, double[] externalU2, bool doNegativeGain, bool doAssertResult=true, double[] yset=null, int pidInputIdx = 0, - double processGainAllowedOffsetPrc=10, bool doAddBadData = false) + double processGainAllowedOffsetPrc=10, bool doAddBadData = false,bool isStatic = true) { var usedProcParameters = trueProcessModel.GetModelParameters().CreateCopy(); var usedProcessModel = new UnitModel(usedProcParameters,"UsedProcessModel"); @@ -332,7 +336,7 @@ public void GenericMISODisturbanceTest (UnitModel trueProcessModel, double[] tr if (doAssertResult) { CommonPlotAndAsserts(pidDataSet, estDisturbance, trueDisturbance, identifiedModel, - trueProcessModel, processGainAllowedOffsetPrc); + trueProcessModel, processGainAllowedOffsetPrc, 30,isStatic); } } diff --git a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs index 188305e..08a8c16 100644 --- a/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs +++ b/TimeSeriesAnalysis.Tests/Tests/ClosedLoopIdentifierTester_StaticSISO.cs @@ -53,7 +53,7 @@ public void StepDisturbanceANDSetpointStep(double distStepAmplitude, double yset var trueDisturbance = TimeSeriesCreator.Step(160, N, 0, distStepAmplitude); var yset = TimeSeriesCreator.Step(50, N, 50, 50+ ysetStepAmplitude);//do step before disturbance CluiCommonTests.GenericDisturbanceTest(new UnitModel(locParameters, "Process"), trueDisturbance, - false, true, yset, precisionPrc); + false, true, yset, precisionPrc,false, true); } /* @@ -62,35 +62,18 @@ I think the reason this test struggles is that closedloopestimator tries to find this may not be as good an assumption as for the step disturbances considered in other tests. */ - [TestCase(5, 1.0), NonParallelizable] - [TestCase(1, 1.0)] - [TestCase(1, 5.0 )]// this only works when the step change is much bigger than the disturbance + [TestCase(5, 1.0,5 )] + [TestCase(1, 1.0,5)] + [TestCase(1, 5.0,5 )]// this only works when the step change is much bigger than the disturbance - public void SinusDisturbanceANDSetpointStep(double distSinusAmplitude, double ysetStepAmplitude) + public void SinusDisturbanceANDSetpointStep(double distSinusAmplitude, double ysetStepAmplitude, double gainPrecisionPrc) { - double precisionPrc = 20; - var trueDisturbance = TimeSeriesCreator.Sinus(distSinusAmplitude,N/8,timeBase_s,N ); var yset = TimeSeriesCreator.Step(N/2, N, 50, 50 + ysetStepAmplitude);//do step before disturbance CluiCommonTests.GenericDisturbanceTest(new UnitModel(modelParameters, "Process"), trueDisturbance, - false, true, yset, precisionPrc,false, isStatic); + false, true, yset, gainPrecisionPrc,false, isStatic); } - // does not work in the static case, no sens in doing dynamic case too - - /* [TestCase(2, Category = "NotWorking_AcceptanceTest")] - [TestCase(1, Category = "NotWorking_AcceptanceTest")] - [TestCase(0.5, Category = "NotWorking_AcceptanceTest")] - public void SinusDisturbance(double distSinusAmplitude) - { - double precisionPrc = 20; - bool doAddBadData = false; - var trueDisturbance = TimeSeriesCreator.Sinus(distSinusAmplitude, N / 8, timeBase_s, N); - var yset = TimeSeriesCreator.Step(N / 2, N, 50, 50 ); - CluiCommonTests.GenericDisturbanceTest(new UnitModel(modelParameters, "Process"), trueDisturbance, - false, true, yset, precisionPrc,doAddBadData,isStatic); - }*/ - // 0.25: saturates the controller // gain of 5 starts giving huge oscillations... @@ -140,15 +123,11 @@ public void RandomWalkDisturbance(double procGain, double distAmplitude, double false, true, yset, precisionPrc,doBadData,isStatic); } - // these tests seem to run well when run individuall, but when "run all" they seem to fail - // this is likely related to a test architecture issue, not to anything wiht the actual algorithm. - - [TestCase(5, 1.0),NonParallelizable]// most difficult - public void StepDisturbanceANDSetpointSinus(double distStepAmplitude, double ysetStepAmplitude) + [TestCase(5, 1.0,5)] + public void StepDisturbanceANDSetpointSinus(double distStepAmplitude, double ysetStepAmplitude, double precisionPrc ) { - Vec vec = new Vec(); - - double precisionPrc = 20;// works when run indivdually + var vec = new Vec(); + var locParameters = new UnitParameters { TimeConstant_s = 0, @@ -160,19 +139,18 @@ public void StepDisturbanceANDSetpointSinus(double distStepAmplitude, double yse var yset = vec.Add(TimeSeriesCreator.Sinus(ysetStepAmplitude, N / 4, timeBase_s, N),TimeSeriesCreator.Constant(50,N)); CluiCommonTests.GenericDisturbanceTest(new UnitModel(locParameters, "StaticProcess"), trueDisturbance, - false, true, yset, precisionPrc); + false, true, yset, precisionPrc,false, true); } - [TestCase(5)] - [TestCase(-5)] - public void LongStepDisturbance_EstimatesOk(double stepAmplitude) + [TestCase(5,5)] + [TestCase(-5,5)] + public void LongStepDisturbance_EstimatesOk(double stepAmplitude, double gainPrecisionPrc) { bool doInvertGain = false; - // Shared.EnablePlots(); N = 1000; var trueDisturbance = TimeSeriesCreator.Step(100, N, 0, stepAmplitude); - CluiCommonTests.GenericDisturbanceTest(new UnitModel(modelParameters, "StaticProcess"), trueDisturbance,doInvertGain); - // Shared.DisablePlots(); + CluiCommonTests.GenericDisturbanceTest(new UnitModel(modelParameters, "StaticProcess"), trueDisturbance,doInvertGain,true, null, gainPrecisionPrc, false, true); + } @@ -210,7 +188,7 @@ public void StepDisturbance_EstimatesOk(double stepAmplitude, double processGain //Shared.EnablePlots(); var trueDisturbance = TimeSeriesCreator.Step(100, N, 0, stepAmplitude); CluiCommonTests.GenericDisturbanceTest(new UnitModel(modelParameters, "Process"), trueDisturbance, - doNegativeGain,true,null, processGainAllowedOffsetPrc); + doNegativeGain,true,null, processGainAllowedOffsetPrc,false, true); //Shared.DisablePlots(); }