------------------- KALMANF.NEW follows --------------------
[Hans Blom, 950522]
Bill, we made a serious error in the interpretation of my demo.
Unmodeled dynamics ARE controlled away by a model-based con-
troller! I have to set some things straight.
I've just had the time to check your statement that a sine wave,
added as "noise" (arbitrary variations) to the "world" will not
be controlled away. I had silently wondered how that could be
true, and I had taken your word for it that it was. It isn't.
Unmodelled dynamics IS resisted, as the following modified demo
shows.
I set ft = 0.1, and I generate a sine wave disturbance with this
amplitude.I also set pff = 0.01, the square of ft. The world is
now more complex than can be modelled.
I set the observation noise vt = 0.0, with pvv = 0.000001 (zero
gives an obscure error message in the plot routines).
The model doesn't know anything about the sine wave, of course.
Yet control eliminates it. Run the demo and see. Verify that I did
not change the code in any other way, except setting the "kludge
factor" delta = 1.5. This forces somewhat faster forgetting. Leav-
ing delta at its old value works all right as well, but generates
more (but smaller) divergences (indicated by the vertical red
lines).
Note that the added sine wave is NOT eliminated while the control-
ler runs blind, as expected.
Can you tell me what you did (wrong) in order to reach your result
and conclusion? I presume -- and my reconstruction demonstrates --
that it was by setting pvv rather than pff to the sine's ampli-
tude. The term pvv tells the model how much OBSERVATION (sensor)
noise is to be expected, i.e. how seriously to take the measure-
ments. Setting pvv to zero says that the measurements are to be
trusted completely. Setting pvv = ft^2 says that the sine wave
disturbance is due to sensor noise and hence must be filtered out
FROM THE MEASUREMENTS because it does not reflect anything in the
"world".
You can verify this by adding a sine wave to y, and setting pvv
equal to its amplitude squared: now the sine wave is disregarded/
filtered out (almost) completely.
This illustrates the difference between ft and vt in a dramatic way.
The term ft is OF the "world", the term vt denotes a faulty PERCEPT-
ION of the world.
Notice the result of the parameter estimates of a, b and c. There
is no convergence anymore to the "real" values of the "world".
Convergence has become impossible because world and model now have
a different structure. Because gamma and delta are set to values
greater than zero, forgetting forces continuous relearning. And
the relearning attempts to subsume the sine wave into the model.
Therefore the parameter estimates become rather meaningless. But
control is fine!
The demo now shows model-based control at its worst: the model is
not complex enough, parameter values do not converge and at all
times have large uncertainties, and keep fluctuating. Yet, thanks
to the forgetting, the model keeps readjusting well enough. The
result is that no certainty or "long term knowledge" develops.
This is an indication that the model is insufficiently complex
and must be expanded. It surprises me that even this limited model
is robust enough to be able to control so well.
I probably introduced the misunderstanding myself when distributing
the code with inappropriate settings for pff and pvv. When doing
new experiments:
set pff to ft^2 (or a little bit more, to make things more stable)
set pvv to vt^2 (or a little bit more, to make things more stable)
Setting pff and pvv to slightly higher values (a factor of 2 to 10)
only makes convergence slower. Setting pff and/or pvv too low makes
for all sorts of problems.
Nobody else discovered this? Does this mean that only Bill Powers
experimented with my demo? Or that my explanation was that unclear?
By the way, "forgetting" in my demo, including the new one, may be
far too drastic and/or done in an inappropriate way. Verify this by
removing "forgetting": just comment out the call to procedure check_
resididual in the main program and see what happens. I still do not
have a good understanding on how to implement "optimal" forgetting.
Here is the demo with the new settings:
program kalman_filter_simulation_with_sinewave_disturbance;
{directions for use of the program:
compile with Turbo Pascal and run
(must have EGAVGA.BGI in current directory!)
press SPACE to view the next batch of runs
press RETURN to end simulation}
uses
crt, graph;
const
minX = 40; {start of plots at this x}
var
run: longint; {increments each sample}
batch: integer; {increments each batch}
t: integer; {step counter within batch}
{system variables}
xt, at, bt, ct: real; {system parameter values}
ft: real; {system noise amplitude}
vt: real; {observation noise amplitude}
xopt: real; {optimal x, setpoint}
{model variables}
x, a, b, c: real; {model parameter values}
pxx, pxa, pxb, pxc,
pax, paa, pab, pac,
pbx, pba, pbb, pbc,
pcx, pca, pcb, pcc: real; {and their covariance matrix}
pff: real; {assumed system noise variance}
pvv: real; {assumed observ. noise variance}
u: real; {action}
y: real; {observation}
r: real; {residual = prediction error}
pyy: real; {expected variance of residual}
qaa, qbb, qcc: real; {expected variation on change}
{screen variables}
MaxX, MaxY: integer; {screen dimensions}
xmin, xmax,
amin, amax,
bmin, bmax,
cmin, cmax,
umin, umax: real; {display scaling factors}
{procedures for plotting}
procedure var_scales;
{define minima and maxima for plotting; change if you like}
begin
xmin := 8.0; xmax := 12.0;
amin := 0.8; amax := 1.0;
bmin := -0.2; bmax := 0.2;
cmin := 0.0; cmax := 2.0;
umin := -6.0; umax := 6.0
end;
procedure init_graphics; {initialize graphics}
var
grDriver,
grMode,
ErrCode,
i, j: integer;
s: string [50];
begin
grDriver := Detect;
InitGraph (grDriver, grMode, '');
ErrCode := GraphResult;
if ErrCode <> grOk then
begin
WriteLn('Graphics error:', GraphErrorMsg (ErrCode));
halt
end;
setcolor (white);
{hope this works on other than EGA/VGA displays!}
maxX := GetMaxX; maxY := GetMaxY - 20;
line (minX, 0, minX, maxY);
line (minX, maxY, maxX, maxY);
line (maxX, maxY, maxX, 0);
for i := 0 to 4 do
begin
line (minX, i * maxY div 5, maxX, i * maxY div 5);
moveto (0, i * maxY div 5 + maxY div 10);
case i of
0: outtext (' x');
1: outtext (' a');
2: outtext (' b');
3: outtext (' c');
4: outtext (' u');
end
end;
for i := 0 to 4 do
begin
moveto (0, i * maxY div 5 + 5);
case i of
0: begin str (xmax:4:1, s); outtext (s); end;
1: begin str (amax:4:1, s); outtext (s); end;
2: begin str (bmax:4:1, s); outtext (s); end;
3: begin str (cmax:4:1, s); outtext (s); end;
4: begin str (umax:4:1, s); outtext (s); end;
end
end;
for i := 1 to 5 do
begin
moveto (0, i * maxY div 5 - 10);
case i of
1: begin str (xmin:4:1, s); outtext (s); end;
2: begin str (amin:4:1, s); outtext (s); end;
3: begin str (bmin:4:1, s); outtext (s); end;
4: begin str (cmin:4:1, s); outtext (s); end;
5: begin str (umin:4:1, s); outtext (s); end;
end
end;
{print info on bottom line}
str (run, s); s := 'run=' + s;
moveto (minX, maxY + 10); outtext (s);
str (batch, s); s := 'batch=' + s + '; SPACE to continue, RETURN to stop';
moveto ((maxX + minX - 8*length (s)) div 2, maxY+10); outtext (s);
str (run+maxX-minX, s); s := 'run=' + s;
moveto (maxX - 8 * length (s), maxY + 10); outtext (s)
end;
procedure plotdot (x, xmin, xmax, w: real; c: integer);
{plot variable x scaled between xmin and xmax in window w, color c}
begin
if x > xmax then x := xmax else if x < xmin then x := xmin;
putpixel (t, round(0.2 * maxY * (w - (x - xmin)/(xmax - xmin))), c)
end;
procedure plot_vars; {plot all interesting variables}
begin
{window 1}
plotdot (x + sqrt (pxx), xmin, xmax, 1, lightgray);
plotdot (x - sqrt (pxx), xmin, xmax, 1, lightgray);
plotdot (x, xmin, xmax, 1, white); {model's x and 1-sigma bounds}
plotdot (xt, xmin, xmax, 1, yellow); {system's x}
plotdot (xopt, xmin, xmax, 1, green); {setpoint}
plotdot (y, xmin, xmax, 1, lightblue);{observation}
{window 2}
plotdot (a + sqrt (paa), amin, amax, 2, lightgray);
plotdot (a - sqrt (paa), amin, amax, 2, lightgray);
plotdot (a, amin, amax, 2, white); {model's a and 1-sigma bounds}
plotdot (at, amin, amax, 2, yellow); {system's a}
{window 3}
plotdot (b + sqrt (pbb), bmin, bmax, 3, lightgray);
plotdot (b - sqrt (pbb), bmin, bmax, 3, lightgray);
plotdot (b, bmin, bmax, 3, white); {model's b and 1-sigma bounds}
plotdot (bt, bmin, bmax, 3, yellow); {system's b}
{window 4}
plotdot (c + sqrt (pcc), cmin, cmax, 4, lightgray);
plotdot (c - sqrt (pcc), cmin, cmax, 4, lightgray);
plotdot (c, cmin, cmax, 4, white); {model's c and 1-sigma bounds}
plotdot (ct, cmin, cmax, 4, yellow); {system's c}
{window 5}
plotdot (u, umin, umax, 5, white) {control}
end;
procedure exit_graphics; {close graphics}
begin
CloseGraph
end;
{auxiliary procedures}
function normal (sd: real): real;
{delivers a zero-average number from a Gaussian distribution}
const
save: real = 0.0;
turn: boolean = false;
var
u, v, c: real;
begin
if turn then {take old number}
normal := sd * save {scale}
else
begin {generate two new numbers}
repeat
u := 2 * random - 1.0;
v := 2 * random - 1.0;
c := u * u + v * v
until c < 1.0; {vector within unit circle}
c := sqrt (- 2.0 * ln (c) / c);
normal := sd * u * c; {scale}
save := v * c
end;
turn := not turn
end;
function chi2test (s: real): boolean;
{test whether s passes chi-square test; returns true if not;
use one of the arrays chi2a, chi2b or chi2c in the code below}
const
chi2a: array [1..30] of real = {90%}
(2.71, 4.61, 6.25, 7.78, 9.24, 10.6, 12.0, 13.4, 14.7, 16.0,
17.3, 18.5, 19.8, 21.1, 22.3, 23.5, 24.8, 26.0, 27.2, 28.4,
29.6, 30.8, 32.0, 33.2, 34.4, 35.6, 36.7, 37.9, 39.1, 40.3);
chi2b: array [1..30] of real = {95%}
(3.84, 5.99, 7.81, 9.49, 11.1, 12.6, 14.1, 15.5, 16.9, 18.3,
19.7, 21.0, 22.4, 23.7, 25.0, 26.3, 27.6, 28.9, 30.1, 31.4,
32.7, 33.9, 35.2, 36.4, 37.7, 38.9, 40.1, 41.3, 42.6, 43.8);
chi2c: array [1..30] of real = {99%}
(6.63, 9.21, 11.3, 13.3, 15.1, 16.8, 18.5, 20.1, 21.7, 23.2,
24.7, 26.2, 27.7, 29.1, 30.6, 32.0, 33.4, 34.8, 36.2, 37.6,
38.9, 40.3, 41.6, 43.0, 44.3, 45.6, 47.0, 48.3, 49.6, 50.9);
shiftr: array [1..30] of real =
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
var
i: byte; {counter}
sum: real; {accumulator}
detect: boolean; {test failed}
begin
{shift new value into shiftr}
for i := 1 to 29 do shiftr [i + 1] := shiftr [i]; shiftr [1] := s;
{accumulate and compare}
sum := 0.0; detect := false;
i := 1;
repeat
sum := sum + shiftr [i];
detect := sum > chi2c [i];
inc (i)
until detect or (i > 30);
chi2test := detect
end;
{extended Kalman filter procedures}
{procedures that describe the behavior of the system}
procedure init_system; {initialize system's parameters}
{these variables describe the system's behavior}
begin
xt := 12.0; {xt := ct + at * xt + bt * u + normal (ft)}
ct := 1.0;
at := 0.9;
bt := 0.1;
ft := 0.1; {system noise amplitude; pff should agree}
vt := 0.0; {observation noise amplitude; pvv should agree}
u := 0.0 {action}
end;
procedure step_system; {execute one step of system}
{instead of normal (ft) try also ft * sin (run/5)}
begin
xt := ct + at * xt + bt * u + ft * sin (run/5) {<----- sine wave}
end;
procedure change_system; {change system's parameters whenever wanted}
begin
case run of
1250: begin ct := 0.8; at := 0.95; bt := 0.15 end;
1500: begin ct := 1.0; at := 0.90; bt := 0.10 end;
2050: begin ct := 0.8; at := 0.95; bt := 0.15 end;
2150: begin ct := 1.0; at := 0.90; bt := 0.10 end;
2600: begin bt := -0.10 end; {sign reversal only}
end
end;
{procedures that describe the observation}
procedure observe; {do one observation}
begin
y := xt + normal (vt)
end;
procedure go_blind; {manipulate pvv whenever wanted}
begin
pvv := 0.000001; {assumed variance of vt}
{pvv is assumed observation noise; 10000 is essentially blind}
if (run > 700) and (run < 900) or
(run > 1850) and (run < 2350) then
begin y := 10; pvv := 1000000 end {freeze y for clarity only}
end;
{procedures that update the model}
procedure init_model; {initialize model variables}
begin
{parameter estimates}
x := 11.5; {initialize these incorrectly}
c := 1.5;
b := 0.15;
pxx := 4.0; {... but not too incorrectly}
pcc := 4.0; {if you do, see what happens!}
paa := 1.0; {correct: initial paa > (at-a)^2, etc.}
pbb := 1.0;
pax := 0.0; pxa := pax;
pbx := 0.0; pxb := pbx;
pcx := 0.0; pxc := pcx;
pab := 0.0; pba := pab;
pac := 0.0; pca := pac;
pbc := 0.0; pcb := pbc;
pff := 0.01; {initial estimate of variance of system noise ft}
qaa := 0.1; {if a changes, by how much?}
qbb := 0.1;
qcc := 0.1;
end;
procedure step_model;
{compute the predicted value of the parameter vector and
its covariance matrix at the next sample interval}
begin
pxx := pcx + a * pxx + x * pax + u * pbx;
pxa := pca + a * pax + x * paa + u * pba;
pxb := pcb + a * pbx + x * pab + u * pbb;
pxc := pcc + a * pcx + x * pac + u * pbc;
pxx := pxc + a * pxx + x * pxa + u * pxb
+ paa * pxx + pax * pax + pff;
x := c + a * x + pax + b * u;
pax := pxa;
pbx := pxb;
pcx := pxc
end;
procedure check_residual; {heuristics for model parameter changes;
this is a collection of heuristics (i.e. procedures not well founded
in a nice theory) that all have the purpose of keeping the model's
estimate of the variance of the residuals consistent with the
actually observed residuals; mainly required in order to be able to
track c.q. readjust to changes of the parameters of the system; in
a noisefree system also to prevent accumulated computational errors}
const
delta = 1.5; {update variances; must be >= 1.0}
gamma = 0.1; {update variances; must be >= 0.0}
begin
{readjust for computational errors}
if pxx < 0 then pxx := - 4.0 * pxx; {pxx etc should not be negative}
if pcc < 0 then pcc := - 4.0 * pcc;
if paa < 0 then paa := - 4.0 * paa;
if pbb < 0 then pbb := - 4.0 * pbb;
{a too large residual should not happen; it indicates an outlier;
increase the modelled uncertainties whenever an outlier is observed;
do this in such a way that consistent occurrences of outliers suf-
ficiently increase the uncertainties, but a single occurrence does
not have too much influence; the following code works fairly well}
r := y - x; {residual = prediction error}
pyy := pxx + pvv; {expected value of residual's variance}
if chi2test (r * r / pyy) then {increase uncertainty}
begin
line (t, 0, t, maxY); {plot outlier position}
pff := delta * pff; {increase uncertainties}
paa := delta * paa + gamma * qaa;
pbb := delta * pbb + gamma * qbb;
pcc := delta * pcc + gamma * qcc;
pxx := delta * pxx;
end
else
pff := 0.99 * pff {decrease estimate of system noise}
end;
procedure process_observation; {process one observation y}
{combine prediction and observation into one probability function}
var
t: real;
begin
t := (y - x) / (pxx + pvv);
x := x + pxx * t;
c := c + pcx * t;
b := b + pbx * t;
t := 1.0 / (pxx + pvv);
pxa := pxa - pxa * pxx * t;
pxb := pxb - pxb * pxx * t;
pxc := pxc - pxc * pxx * t;
pxx := pxx - pxx * pxx * t;
paa := paa - pax * pax * t;
pab := pab - pax * pbx * t;
pac := pac - pax * pcx * t;
pbb := pbb - pbx * pbx * t;
pbc := pbc - pbx * pcx * t;
pcc := pcc - pcx * pcx * t;
pax := pxa;
pbx := pxb;
pba := pab;
pcx := pxc;
pca := pac;
pcb := pbc;
end;
{procedures that define how to control}
procedure setxopt; {define setpoint}
begin
xopt := 10.0 + sin (run / 25.0)
end;
procedure control; {compute control; could be improved considerably;
{note that the computation of the control u is based entirely upon
the model's parameters only; NOT on those of the system; NOT on the
observation y}
begin
if b <> 0.0 then
u := (xopt - c - a * x - pax) / b
else
u := 0.0 {if b = 0 then no control is possible, so don't}
end;
{main program}
begin
init_system; {initialize system's parameters}
init_model; {initialize model's parameters}
var_scales; {set bounds for plotting in windows}
run := 0; batch := 0; {startup}
repeat
inc (batch); {next batch of runs}
t := minX - 1; {t = screen x-position for plotting}
init_graphics; {setup axes etc. on screen}
setcolor (red); {outliers signalled by red vertical line}
repeat
inc (run); inc (t); {next run}
step_model; {predict across one sample}
step_system; {... the real thing also}
change_system; {system parameter changes ?}
observe; {observe the system}
go_blind; {missing observation ?}
check_residual; {observation fits prediction ?}
process_observation; {collapse of the wave function}
setxopt; {(re)define setpoint}
control; {compute control}
plot_vars; {plot interesting variables}
until t = maxX; {display screen full ?}
repeat until keypressed;
until readkey = chr (13); {stop when RETURN pressed}
exit_graphics {shut down graphics system}
end.
···
a := 0.8;
a := a + pax * t;