Time the variable in dynamic optimization using orthogonal collocation is supposed to act

Question

I am having difficulties wrapping my head around when exactly does the manipulated variable act in orthogonal collocation.

I have the following problem:

$$ \min x_2(tf) \\ \frac{dx_1}{dt}=u\\\frac{dx_2}{dt}=x_1^2 + u^2 \\ x_0 = [1, 0]^T\\t_f = 1$$

I have used three collocation points and programmed it in Pyomo. The code is:

T = [0, 0.155051, 0.644949, 1]
t1 = T[1]
t2 = T[2]
t3 = (1)
M1 = np.array([[1, 2*t1, 3*t1**2],
               [1, 2*t2,  3*t2**2],
               [1, 2*t3, 3*t3**2]])
M2 = np.array([[t1, t1**2, t1**3],
               [t2, t2**2, t2**3],
               [t3, t3**2, t3**3]]) 
OMEGA = np.linalg.inv(M1@(np.linalg.inv(M2)))

model = pyo.ConcreteModel()
model.i = pyo.Set(initialize = [1, 2, 3])
model.k = pyo.Set(initialize = [1, 2, 3])
model.j = pyo.SetOf(model.k)

model.x1 = pyo.Var(model.i, model.k)
model.dx1 = pyo.Var(model.i, model.k)
model.x2 = pyo.Var(model.i, model.k)
model.dx2 = pyo.Var(model.i, model.k)
model.x10 = pyo.Var(model.i)
model.x20 = pyo.Var(model.i)
model.u = pyo.Var(model.i, model.k)
OMEGA_dict = {}
for n,j in enumerate(model.j):
    for nn,k in enumerate(model.k):
        OMEGA_dict[j,k] = OMEGA[n,nn]
model.OMEGA = pyo.Param(model.j, model.k, initialize = OMEGA_dict)

def state_discretizedx1(m, i, k):
    return m.x1[i,k] == m.x10[i] + 0.33333*sum(model.OMEGA[k,j]*m.dx1[i,j] for j in m.k)
model.eq1 = pyo.Constraint(model.i, model.k, rule = state_discretizedx1)

def state_discretizedx2(m, i, k):
    return m.x2[i,k] == m.x20[i] + 0.33333*sum(model.OMEGA[k,j]*m.dx2[i,j] for j in m.k)
model.eq2 = pyo.Constraint(model.i, model.k, rule = state_discretizedx2)

def dx1dt(m, i, k):
    return m.dx1[i,k] == m.u[i,k]
def dx2dt(m,i,k):
    return m.dx2[i,k] == m.x1[i,k]**2 + m.u[i,k]**2
model.eq3 = pyo.Constraint(model.i, model.k, rule=dx1dt)
model.eq4 = pyo.Constraint(model.i, model.k, rule=dx2dt)

def continuityx1(m, i, k):
    if i > 1:
        return m.x10[i] == m.x1[i-1, 3]
    else:
        return pyo.Constraint.Skip
model.eq5 = pyo.Constraint(model.i, model.k, rule = continuityx1) 

def continuityx2(m, i, k):
    if i > 1:
        return m.x20[i] == m.x2[i-1, 3]
    else:
        return pyo.Constraint.Skip
model.eq6 = pyo.Constraint(model.i, model.k, rule = continuityx2) 


model.obj = pyo.Objective(expr = model.x2[len(model.i),3])

model.x10[1].fix(1)
model.x20[1].fix(0)

It solves perfectly, but I have problems interpreting the result. For example, the manipulated variable $u$, only has values for the collocation points $1,2,3$. Therefore, for example, the solution says that $u_{i1,k1} = -0.711$.

However, I am not sure when this action should be taken. $i1,k1$ means at time 0.155, according to the Radau roots used. What would be the action in the period $[0, 0.155)$ then? How does this work if I wanted to actually apply the solution?

Thank you