Skip to content

Commit

Permalink
- CLDI wip and improving documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
Steinar Elgsæter committed Dec 9, 2024
1 parent 9a2388e commit 35f039a
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 16 deletions.
42 changes: 32 additions & 10 deletions Dynamic/Identification/ClosedLoopUnitIdentifier.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
////////////////////////
Expand All @@ -54,6 +54,8 @@ public class ClosedLoopUnitIdentifier
/// <returns>The unit model, with the name of the newly created disturbance added to the additiveInputSignals</returns>
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)
Expand Down Expand Up @@ -125,15 +127,21 @@ 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;
idDisturbancesList.Add(distIdResult1);
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];
Expand All @@ -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
Expand Down Expand Up @@ -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));
}
}

Expand Down
11 changes: 10 additions & 1 deletion Dynamic/Identification/DisturbanceIdentifier.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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());
Expand Down Expand Up @@ -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;
Expand All @@ -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<double>.Fill(1, nGains);
Expand Down Expand Up @@ -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
Expand Down
23 changes: 23 additions & 0 deletions TimeSeriesAnalysis/Vec.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -721,6 +721,29 @@ public double[] Multiply(double[] array1, double[] array2)
return retVal;
}

/// <summary>
/// Returns the mean value of array1. while ignoring given indices
/// </summary>
public double? Mean(double[] array1, List<int> 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;
}






Expand Down
29 changes: 24 additions & 5 deletions docs/articles/sysid_disturbance.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -114,18 +114,37 @@ a linear static model essentially boils down to estimating the process gain.

![init](./images/sysid_disturbance_init.png)

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

Expand Down

0 comments on commit 35f039a

Please sign in to comment.