diff --git a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs
index 2ab0760..3304252 100644
--- a/Dynamic/Identification/ClosedLoopUnitIdentifier.cs
+++ b/Dynamic/Identification/ClosedLoopUnitIdentifier.cs
@@ -37,7 +37,7 @@ public class ClosedLoopUnitIdentifier
// NB!! These three are somewhat "magic numbers", that need to be changed only after
// testing over a wide array of cases
const int firstPassNumIterations = 30;//TODO: change back to 50
- const int secondPassNumIterations = 20;
+ const int secondPassNumIterations = 0;
const double initalGuessFactor_higherbound = 2.5;// 2 is a bit low, should be a bit higher
const int nDigits = 5; //number of significant digits in results.
////////////////////////
@@ -54,6 +54,8 @@ public class ClosedLoopUnitIdentifier
/// The unit model, with the name of the newly created disturbance added to the additiveInputSignals
public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters pidParams = null, int pidInputIdx = 0)
{
+ bool doConsoleDebugOut = true;
+
bool wasGainGlobalSearchDone = false;
bool doTimeDelayEstOnRun1 = false;
if (dataSet.Y_setpoint == null || dataSet.Y_meas == null || dataSet.U == null)
@@ -125,6 +127,9 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters
var distIdResult1 = DisturbanceIdentifier.EstimateDisturbance
(dataSetRun1, null, pidInputIdx, pidParams);
+ if (doConsoleDebugOut)
+ Console.WriteLine("Run1,step1: " + distIdResult1.estPidProcessGain.ToString("F3", CultureInfo.InvariantCulture));
+
dataSetRun1.D = distIdResult1.d_est;
var unitModel_run1 = UnitIdentifier.IdentifyLinearAndStatic(ref dataSetRun1, fittingSpecs, doTimeDelayEstOnRun1);
unitModel_run1.modelParameters.LinearGainUnc = null;
@@ -132,8 +137,11 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters
idUnitModelsList.Add(unitModel_run1);
isOK = ClosedLoopSim(dataSetRun1, unitModel_run1.GetModelParameters(), pidParams, distIdResult1.d_est, "run1");
+ if (doConsoleDebugOut)
+ Console.WriteLine("Run1,ident: " + unitModel_run1.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture));
+
// run1, "step2" : "global search" for linear pid-gainsgains
- if(unitModel_run1.modelParameters.GetProcessGains()!= null)
+ if (unitModel_run1.modelParameters.GetProcessGains()!= null)
{
double pidProcessInputInitalGainEstimate = unitModel_run1.modelParameters.GetProcessGains()[pidInputIdx];
@@ -150,27 +158,39 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters
var max_gain = Math.Abs(pidProcessInputInitalGainEstimate * initalGuessFactor_higherbound);
var min_gain = - max_gain;
- // min_gain = 0; // when debugging, it might be advantageous to set min_gain equal to the known true value
+ if (doConsoleDebugOut)
+ {
+ Console.WriteLine("Run1,GS Higher bound: "+ max_gain.ToString("F3", CultureInfo.InvariantCulture) );
+ }
+
+ // min_gain = 0; // when debugging, it might be advantageous to set min_gain equal to the known true value
// first pass(wider grid with larger grid size)
- var retPass1 = GlobalSearchLinearPidGain(dataSet, pidParams, pidInputIdx,
+ var retGlobalSearch1 = GlobalSearchLinearPidGain(dataSet, pidParams, pidInputIdx,
unitModel_run1, pidProcessInputInitalGainEstimate,
min_gain, max_gain, fittingSpecs,firstPassNumIterations);
- var bestUnitModel = retPass1.Item1;
+ var bestUnitModel = retGlobalSearch1.Item1;
+
+ if (doConsoleDebugOut && retGlobalSearch1.Item1 != null)
+ Console.WriteLine("Run1,GS1: " + retGlobalSearch1.Item1.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture));
+ else
+ Console.WriteLine("Run1,GS1: FAILED");
if (bestUnitModel != null)
{
// second pass(finer grid around best result of first pass)
- if (retPass1.Item1.modelParameters.Fitting.WasAbleToIdentify && secondPassNumIterations > 0)
+ if (retGlobalSearch1.Item1.modelParameters.Fitting.WasAbleToIdentify && secondPassNumIterations > 0)
{
const int WIDTH_OF_SEARCH_PASS2 = 3;
- var gainPass1 = retPass1.Item1.modelParameters.LinearGains[pidInputIdx];
- var retPass2 = GlobalSearchLinearPidGain(dataSet, pidParams, pidInputIdx,
- retPass1.Item1, gainPass1, gainPass1 - retPass1.Item2* WIDTH_OF_SEARCH_PASS2, gainPass1 + retPass1.Item2* WIDTH_OF_SEARCH_PASS2,
+ 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,
fittingSpecs,secondPassNumIterations);
- bestUnitModel = retPass2.Item1;
+ bestUnitModel = retGlobalSearch2.Item1;
wasGainGlobalSearchDone = true;
+ if (doConsoleDebugOut)
+ Console.WriteLine("Run1,GS2: " + retGlobalSearch2.Item1.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture));
}
}
// add the "best" model to be used in the next model run
@@ -324,6 +344,8 @@ public static (UnitModel, double[]) Identify(UnitDataSet dataSet, PidParameters
var distIdResult_step4 = DisturbanceIdentifier.EstimateDisturbance
(dataSetRun2, step2Model, pidInputIdx, pidParams);
idDisturbancesList.Add(distIdResult_step4);
+ if (doConsoleDebugOut)
+ Console.WriteLine("Run2 " + step2Model.GetModelParameters().LinearGains.First().ToString("F3", CultureInfo.InvariantCulture));
}
}
diff --git a/Dynamic/Identification/DisturbanceIdentifier.cs b/Dynamic/Identification/DisturbanceIdentifier.cs
index 338f522..3eb7330 100644
--- a/Dynamic/Identification/DisturbanceIdentifier.cs
+++ b/Dynamic/Identification/DisturbanceIdentifier.cs
@@ -256,7 +256,6 @@ private static UnitModel EstimateClosedLoopProcessGain(UnitDataSet unitDataSet,
var unitModel = new UnitModel();
var vec = new Vec(unitDataSet.BadDataID);
- //var result = new DisturbanceIdResult(unitDataSet);
if (unitDataSet.Y_setpoint == null || unitDataSet.Y_meas == null || unitDataSet.U == null)
{
return null;
@@ -277,6 +276,7 @@ private static UnitModel EstimateClosedLoopProcessGain(UnitDataSet unitDataSet,
// but works better than candiate 2 when disturbance is a step
double FilterTc_s = 0;
+
// initalizaing(rough estimate):
{
LowPass lowPass = new LowPass(unitDataSet.GetTimeBase());
@@ -316,6 +316,12 @@ private static UnitModel EstimateClosedLoopProcessGain(UnitDataSet unitDataSet,
double maxU = vec.Max(vec.Abs(uFiltered), unitDataSet.IndicesToIgnore); // sensitive to output noise/controller overshoot
double minU = vec.Min(vec.Abs(uFiltered), unitDataSet.IndicesToIgnore); // sensitive to output noise/controller overshoot
estPidInputProcessGain = pidInput_processGainSign * maxDE / (maxU - minU);
+
+ // double avgDE = vec.Mean(vec.Abs(eFiltered), unitDataSet.IndicesToIgnore).Value;
+ // double avgU = vec.Mean(vec.Abs(pidInput_deltaU), unitDataSet.IndicesToIgnore).Value;
+ // double estPidInputProcessGainUB = avgDE / avgU;
+ // Console.WriteLine("process gain upper og lower boundbound:"+estPidInputProcessGainUB+ " uncertainty:" + Math.Abs(estPidInputProcessGain- estPidInputProcessGainUB));
+
}
int indexOfFirstGoodValue = 0;
@@ -337,6 +343,8 @@ private static UnitModel EstimateClosedLoopProcessGain(UnitDataSet unitDataSet,
var unitParamters = new UnitParameters();
unitParamters.LinearGains = new double[nGains];
unitParamters.LinearGains[pidInputIdx] = estPidInputProcessGain;
+ // unitParamters.LinearGainUnc[pidInputIdx] = estPidInputProcessGain-est;
+
unitParamters.U0 = new double[nGains];
unitParamters.U0[pidInputIdx] = pidInput_u0[pidInputIdx];
unitParamters.UNorm = Vec.Fill(1, nGains);
@@ -432,6 +440,7 @@ public static DisturbanceIdResult EstimateDisturbance(UnitDataSet unitDataSet,
result.d_est = d_est;
result.d_LF = d_LF;
result.d_HF = d_HF;
+ result.estPidProcessGain = unitModel.GetModelParameters().LinearGains.ElementAt(pidInputIdx);
result.adjustedUnitDataSet = unitDataSet_adjusted;
// END STEP 2
diff --git a/TimeSeriesAnalysis/Vec.cs b/TimeSeriesAnalysis/Vec.cs
index afc3f6f..2ea5a3f 100644
--- a/TimeSeriesAnalysis/Vec.cs
+++ b/TimeSeriesAnalysis/Vec.cs
@@ -721,6 +721,29 @@ public double[] Multiply(double[] array1, double[] array2)
return retVal;
}
+ ///
+ /// Returns the mean value of array1. while ignoring given indices
+ ///
+ public double? Mean(double[] array1, List indToIgnore = null)
+ {
+ if (array1 == null)
+ return null;
+ double retVal = 0;
+ double N = 0;
+ var indices = Index.InverseIndices(array1.Length, indToIgnore);
+ foreach (var ind in indices)
+ {
+ if (!IsNaN(array1[ind]))
+ {
+ N += 1;
+ retVal = retVal * (N - 1) / N + array1[ind] * 1 / N;
+ }
+ }
+ return retVal;
+ }
+
+
+
diff --git a/docs/articles/sysid_disturbance.md b/docs/articles/sysid_disturbance.md
index 8ebe0d8..c3026c3 100644
--- a/docs/articles/sysid_disturbance.md
+++ b/docs/articles/sysid_disturbance.md
@@ -105,7 +105,7 @@ to the algorithm, so that the algorith can infer about the control error ``e``.
> is also true, what appear to be negative process gains at first sight may in
> closed-loop be positive process gains.
-### First, model-free estimate
+### First, model-free estimate of the process gain
A model-free estimate of the disturbance is required to initialize
subsequent sequential estimation.
@@ -114,18 +114,37 @@ a linear static model essentially boils down to estimating the process gain.
-This first estimate of the process gain ``G`` in a linear model ``y = G x u``
+This first estimate of the process gain ``G_0`` in a linear model ``y = G_0 x u``
is found by the approximation
```
-G = max(e)/(max(u)-min(u))
+G_0 = max(e)/(max(u)-min(u))
```
+
+The PID-controller usually has an integral effect which causes it to be somewhat delayed in
+responding to a disturbance, and the process also can have time delay and/or a time constant
+that means that the deviation ``e`` is brought back to zero over a certain period of time.
+The idea of creating an inital estimate withh ``min`` and ``max`` values is that it circumvents the lack
+of knowledge of the dynamics at this early stage of estimaton.
+
+It has been observed in unit tests that this estiamte in some cases is spot on the actual gain, such as when
+the disturbance is a perfect step.
+
+Remember, it is not needed for the inital heuritic guess of the process gain to be completely accurate, because
+subsequent steps in the search algorithm can try to improve on this estiamte.
+
+Once an inital estiamte of the process gain has been made, it is possible to estimate the disturbance vector, which is used to initate the below steps.
+
### Pid gain global search
-TODO
+After the inital estiamte of the
+
+
+### Estimating upper-and lower boundds on the process gain
+
### Time constant search algorithm
-TODO
+
### Separating the state of the system from the output