Solving recursive matrix system not fully correct

Question

I'm having a problem with a personal project in which I want to determine the dynamics of a population.

Suppose the following transition matrix is given by

\begin{equation} T = \begin{pmatrix} 0.9 & 0 & 0 \\ 0.05 & 0.95 & 0 \\ 0.05 & 0.05 & 1 \end{pmatrix} \end{equation} The population is normally then determined by \begin{equation} P_{n} = TP_{n-1} \end{equation} However, in this problem the population at $n$ is determined by a different equation involving a sort of manual extra input matrix, i.e. \begin{equation} P_{n} = TP_{n-1} + S \end{equation} Now, this problem will differ from regular transition matrix problems (for as far as I know). After $n$ years, I want to have a pre-determined optimal population, say at $n=f$ (final year). Also, the starting population is of course known. We can now simulate this matrix system up to $n=f$ such that we can see whether we have reached the optimal population, say $P_f$. However, as this is interesing too, I want to solve this recursive matrix system for $S$ such that I find the suitable $S$ matrix such that after $n=f$ years the optimal population is obtained. I tried to do the following

\begin{align*} P_{n} &= AP_{n-1} + S\\ &= A(AP_{n-2} + S) + S \\ &= A^2P_{n-2} + AS + S \\ &= A^2(AP_{n-3} + S) + AS + S \\ &= A^3P_{n-3} + A^2S + AS + S \end{align*} or phrased differently \begin{align*} P_{n} = A^n P_0 + \sum_{j=1}^n A^{j-1}S \end{align*} Then, simply solving for $S$ would yield

\begin{equation} S = \bigg(\sum_{j=1}^n A^{j-1}\bigg)^{-1} \bigg(P_n - A^n P_0\bigg) \end{equation}

However, when I try to simulate my answer, it does not seem to work completely. Suppose I have the following data

\begin{equation} P_0 = \begin{pmatrix} 85 \\ 110 \\ 0 \end{pmatrix} \end{equation}

\begin{equation} T = \begin{pmatrix} 0.9 & 0 & 0 \\ 0.05 & 0.95 & 0 \\ 0.05 & 0.05 & 1 \end{pmatrix} \end{equation}

\begin{equation} P_f = P_4 = \begin{pmatrix} 80 \\ 130 \\ 0 \end{pmatrix} \end{equation}

substituting the latter into our found solution yields

\begin{equation} S = \begin{pmatrix} 7 \\ 10.8 \\ -1 \end{pmatrix} \end{equation} and running a simulation will yield \begin{equation} P_4 = \begin{pmatrix} 80.000 \\ 145.564 \\ 3.68 \end{pmatrix} \neq P_f \end{equation} Only my the entry seems correct. The second does not. When I transpose the matrix, however, the first entry will become incorrect whereas the second entry will become correct. How is this possible? Is there a mistake or is my request simply not possible?

EDIT: My problem has been solved, and a numerical method in python can be found as the answer to this question. However, I also wanted to post my R-code for the people who are interested.

T_matrix = t(matrix(c(0.9,  0,    0,
                      0.05, 0.95, 0,
                      0.05, 0.05, 1), nrow=3))

P = matrix(c(85,110,0), nrow = 3)
GP = matrix(c(80,130,0), nrow = 3)

years <- 4
T_mat <- NULL
for (j in 1:years){
        T_mat[[j]] <- T_matrix%^%(j-1)
}

S  <- (Reduce("+", T_mat) %>% inv()) %*% (GP - T_matrix%^%years %*% P)

pop <- NULL
for(i in 1: (years+1)){
   if (i == 1){
        pop[[1]] <- P
   }else{
        pop[[i]] <- T_matrix %*% pop[[i-1]] + S  
   }
}

final_population <- pop[[years + 1]] %>% print()
```

Answer 1

You must be making some sort of computational error. I get a completely different value for $S$ that does produce the desired value when the iteration is run.

Here is my python script:

import numpy as np

T= np.array([[.9,0,0],[.05,.95,0],[.05,.05,1]])
P0 = np.array([85,110,0])
P4 = np.array([80,130,0])
A = sum(np.linalg.matrix_power(T, j) for j in range(4))
A1 = np.linalg.inv(A)
S=A1@(P4-np.linalg.matrix_power(T,4)@P0)
print("S=",S)

R = P0
for n in range(4):
    R= T@R+S
print("PF=",R)

This produces the output:

S= [7.04608898   6.74717393 -10.04326291]
PF= [8.00000000e+01  1.30000000e+02 -5.32907052e-15]

We see that PF agrees with the given value P4 except for roundoff error.

It's impossible to say what your mistake is without seeing the details of your calculation. The matrix algebra is correct.