Dual Control Example

Here is Matlab code for a simple dual control example. There is no measurement noise, so the x(i) are known perfectly. The control gain b is not known perfectly but with variance var_b(i). There is process noise with variance var_w = 1. The first section below is the top level code, and the second file is a subroutine that computes the cost of an initial action.

Matlab code to plot cost of various initial actions on a 2 step problem:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% x(k+1) = x(k) + b(k)*u(k) + w
% b(k+1) = b(k)
% w is zero mean noise with variance var_w = 1

% x0 is first state
x0 = -1;

% nominal estimate of b remains the same throughout
b = 1;

% initial b parameter variance
var_b0 = 1;

% var_w: variance of process noise w
var_w = 1

% u0 is first action

N = 100;
u0s(N) = 0;
cs(N) = 0;

for i = 1:N
 u0 = (i-1)/(N-1);
 c = cost( u0, x0, b, var_b0, var_w );
 u0s(i) = u0;
 cs(i) = c;
end

plot( u0s, cs );
xlabel('u0')
ylabel('cost')
title('x0=-1,b=1,var_w=1')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Matlab subroutine to compute cost of initial action

function c = cost( u0, x0, b, var_b0, var_w )
% c = cost( u0, x0, b, var_b0, var_w ): compute cost of initially taking
% action u0
% x0 is the starting state,
% b is the prior on the value of the b parameter,
% var_b0 is the initial variance of b
% x(k+1) = x(k) + b(k)*u(k) + w
% b(k+1) = b(k)
% w is zero mean noise with variance var_w

% step indices are appended to variable names:
% x1 = x(1), x2 = x(2), var_b1 = var_b(1), u1 = u(1), ...

% variance of b after first step
var_b1 = var_b0*var_w/(u0*u0*var_b0 + var_w);

% expectation of x1
x1 = x0 + b*u0;

% use cautious control on last step since further learning is not useful
% since no further decisions will be made
u1 = -b*x1/(var_b1 + b*b);

% expectation of x2
x2 = x1 + b*u1;

% Var(x1)
var_x1 = var_b0*u0*u0 + var_w;

% Var(x2)
var_x2 = var_b1*u1*u1 + var_w;

% E(x1^2)
ex1_2 = var_x1 + x1*x1;

% E(x2^2)
ex2_2 = var_x2 + x2*x2;

% cost which is sum of E(x1^2 + x2^2)
c = ex1_2 + ex2_2;

In the figure, notice that the optimum is around 0.6. If this had been just a one step problem, the optimum u0 would have been 0.5 (using cautious control). Having more than one step makes it advantagous to use a bigger u0 to reduce the uncertainty about the parameter b.

Matlab code to plot optimal u versus other variables:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear uarray

x0 = -1;
b = 1;
var_b0 = 1;
var_w = 1;

i = 1;

for x0 = 0.0:-0.1:-2.0
 j = 1;
 for var_b0 = 0.0: 0.1: 5.0
 u = findbestu( x0, b, var_b0, var_w );
 uarray(i,j) = u;
 j = j + 1;
 end
i = i + 1;
end

figure(1)
y = 0.0:-0.1:-2.0;
x = 0.0: 0.1: 5.0;
surf(x,y,uarray)
view(-37.5+180,30)
title('Optimal u0 when varying x0 and Var_b0')
ylabel('x0')
xlabel('Var(b)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Matlab subroutine to find best u

function u = findbestu( x0, b, var_b0, var_w )
% x0 is first state
% b: nominal estimate of b remains the same throughout
% var_b0: initial b parameter variance
% var_w: variance of process noise

% x(k+1) = x(k) + b(k)*u(k) + w
% b(k+1) = b(k)
% v is zero mean noise with variance var_w

% u0 is first action
% c = cost( u0, x0, b, var_b0 )

N = 1001;
u0s(N) = 0;
cs(N) = 0;

u_best = 0;
c_best = 1e30;

for i = 1:N
 u0 = 10*(i-1)/(N-1);
 c = cost( u0, x0, b, var_b0, var_w );
 if ( c < c_best )
  u_best = u0;
  c_best = c;
 end
 u0s(i) = u0;
 cs(i) = c;
end

u = u_best;

Matlab code to plot u0 vs b and var_b0:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear uarray

x0 = -1;
b = 1;
var_b0 = 1;
var_w = 1;

i = 1;

for b = 0.1:0.1:2.0
 j = 1;
 for var_b0 = 0: 0.1: 5.0
  u = findbestu( x0, b, var_b0, var_w );
  uarray(i,j) = u;
  j = j + 1;
 end
i = i + 1;
end

figure(2)
y = 0.1:0.1:2.0;
x = 0: 0.1: 5.0;
surf(x,y,uarray)
title('Optimal u0 when varying b and Var(b)')
ylabel('b')
xlabel('Var(b)')
view(-37.5+180,30)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Matlab code to plot u0 vs var_w and var_b0:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear uarray

x0 = -1;
b = 1;
var_b0 = 1;
var_w = 1;

i = 1;

for var_w = 0.1:0.1:2.5
 j = 1;
 for var_b0 = 0: 0.1: 5.0
  u = findbestu( x0, b, var_b0, var_w );
  uarray(i,j) = u;
  j = j + 1;
 end
i = i + 1;
end

figure(3)
y = 0.1:0.1:2.5;
x = 0: 0.1: 5.0;
surf(x,y,uarray)
title('Optimal u0 when varying var_w and Var(b)')
ylabel('var_w')
xlabel('Var(b)')
view(-37.5+180,30)

% var_w has no effect

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%