Model prediction gives stochastic inputs - Predictive Control

I have a discrete state space model:

$$x(k+1) = Ax(k) + Bu(k)$$ $$y(k) = Cx(k)$$

And I'm trying to compute the predicted inputs. The first thing I do is that I fist create the extended observability matrix $$\Phi$$

$$\Phi = \begin{bmatrix} CA\\ CA^2\\ CA^3\\ \vdots \\ CA^{n-1} \end{bmatrix}$$

Then create the lower traingular topelitz matrix $$\Gamma = \begin{bmatrix} CB & 0 & 0 & 0 &0 \\ CAB & CB & 0 & 0 &0 \\ CA^2B & CAB & CB & 0 & 0\\ \vdots & \vdots &\vdots & \ddots & 0 \\ CA^{n-2} B& CA^{n-3} B & CA^{n-4} B & \dots & CA^{n-j} B \end{bmatrix}$$

Where $$j=n$$

For SISO-case, as in this case. I compute the predicted inputs signals by using:

$$U = (\Gamma)^{-1}(R-\Phi x)$$

For MIMO case. I need to use this equation:

$$U = (\Gamma^T \Gamma)^{-1}\Gamma^T(R-\Phi x)$$

Where $$R$$ is my reference vector and $$x$$ is the initial state. I have matlab code for that.

function [U] = predict (A, B, C, x, N, r)

## Find matrix
PHI = phiMat(A, C, N);
GAMMA = gammaMat(A, B, C, N);
U = inv(GAMMA'*GAMMA)*GAMMA'*(repmat(r, N, 1) -PHI*x);
endfunction

function PHI = phiMat(A, C, N)

## Create the special Observabillity matrix
PHI = [];
for i = 1:N
PHI = vertcat(PHI, C*A^i);
endfor

endfunction

function GAMMA = gammaMat(A, B, C, N)

## Create the lower triangular toeplitz matrix
GAMMA = [];
for i = 1:N
GAMMA = horzcat(GAMMA, vertcat(zeros((i-1)*size(C*A*B, 1), size(C*A*B, 2)),cabMat(A, B, C, N-i+1)));
endfor

endfunction

function CAB = cabMat(A, B, C, N)

## Create the column for the GAMMA matrix
CAB = [];
for i = 0:N-1
CAB = vertcat(CAB, C*A^i*B);
endfor

endfunction


If I have a state space model, discrete with sample time $$h=0.5$$.

A =

0.89559   0.37735
-0.37735   0.51825

B =

0.10441
0.37735

C =

1   0


When I run this code above. I get this result if I choose $$N = 40$$

   239.4510
-301.5661
301.1320
-208.4868
222.4277
-141.9374
166.1560
-94.3562
125.9231
-60.3368
97.1576
-36.0137
76.5909
-18.6233
61.8862
-6.1896
51.3728
2.7002
43.8559
9.0562
38.4815
13.6006
34.6389
16.8497
31.8916
19.1727
29.9273
20.8336
28.5229
22.0211
27.5188
22.8702
26.8009
23.4772
26.2876
23.9113
25.9206
24.2216
25.6582
24.4434


The correct input signals that tracks the reference is 25. So the algorithm computes correct at last. But when I heard prediction, I was told that I should select the first value of the input signals, not the last one. The first value is 239.4510 and the last value is 24.4434.

Questions:

• I'm I wrong according to prediction that I should use a lower triangular toeplitz matrix?

• Should prediction, by solving linear least squares, looks like this? In this case, I think I need to select the last value.

• If I have a model with a pole equals 1, e.g integration. How would this affect this prediction?

For Model Predictive Control by using a QP-solver, then I can use the first value. But it this case...no.

EDIT:

Here is a simulation when we say that the input signals cannot be $$u < 0$$

I have not calculated the formulation by myself. But, I can certainly tell you that you should use the head of $$\boldsymbol{U}$$ instead of its tail.

Using the tail of $$\boldsymbol{U}$$ is definitely wrong. Look at the trend of the optimal control inputs. They start shrinking. The system looks stable. But, if you read the inputs upside-down then the system's input is growing (with oscillation). That's scary and looks unstable.

Ignore the the number 25. Simulate the system as it is and then see how the system behaves.

• Head = first number and tail = last number ?
– DanM
Jul 11, 2019 at 14:46
• @DanielMårtensson, exactly. Jul 12, 2019 at 10:04
• Bur if I use the first number, the dynamical system is going to looks that it's has shaked a lot. If I use the last number, then the system is going to following the reference very smooth.
– DanM
Jul 13, 2019 at 9:07
• @DanielMårtensson, You cannot judge a signal without a simulation. A high oscillation at the end of the prediction horizon is a very bad sign. So, that theory is not going to solve anything. If you want to reduce oscillation at the beginning of prediction horizon, penalize $\Delta u$. Example at this IEEE conf paper. Jul 15, 2019 at 11:23
• You men regularization? Can you explain how? :)
– DanM
Jul 15, 2019 at 12:03