How can we control the output temperature of a multi-input single-output pipe air-flow system?

There might be so many methods to perform that. But, the premium method is to use a nonlinear Model Predictive Control (MPC).

This is a multi-input/multi-output control method which predicts the future based on the available model. Then it optimizes a series of control inputs to minimize a forecast cost function with respect to the constraints.

A nonlinear MPC is like this

begin{cases} J(Delta boldsymbol U|t)=sum_{k=1}^{N_p} ||T_4(t+kT_s|t)-T_4^{text{ref}}||^2+sum_{k=0}^{N_c-1}||Deltaboldsymbol U||^2 Delta boldsymbol U^{star}(t) =text{argmin}_{Delta boldsymbol U} ~~J(Delta boldsymbol U|t) text{Optional constraints} end{cases}

where $x(t'|t)$ means $x$ predicted of time $t'$ at the current time of $t$.

The first term in the cost function is penalizing the error. The second term minimizes the control actions. The constraints could be maximum or minimum temperatures. $N_p$ is called the prediction horizon and $N_c$ is the control horizon and $T_s$ is the sample time.

$boldsymbol U$ is the series of control actions such as

$$ boldsymbol U=begin{bmatrix} T_1(t) T_2(t) T_3(t) dots T_1(t+(N_c-1)T_s) T_2(t+(N_c-1)T_s) T_3(t+(N_c-1)T_s) end{bmatrix} $$

$Delta boldsymbol U$ is the difference of control input $$Delta boldsymbol U(t)=boldsymbol U(t)-boldsymbol U(t-T_s)$$

The control input is updated at each sample time.

If there is no constraint and if the problem is convex the problem becomes far easier to solve. You can use gradient descent method. If it is even linear, it will reduce to Linear Quadratic Regulator (LQR) problem.

For more information about MPC, have a look at my papers in Researchgate.

BTW, MATLAB has MPC toolbox as well as system identification. You can easily use them.

Let us take a lumped mass approach to model the thermal dynamics of the region across two zone, including the convection and conduction and approximate the state $T_4 approx frac{T_2 + T_3}{2}$. Further, assuming that the heat flow is dominated by convection, we can apply the Newton's law of cooling (or heating) to obtain the heat flow rates and therefore the lumped-mass zone temperature dynamics are $dot{T_{i}} = {k}_{i} (T_i - frac{T_{iup}+ T_{idown}}{2})$ where $0 < {k}_{i} = const.$ and the subscripts $up$, $down$ denote the upstream and downstream temperatures. Then, assuming that $k := k_i = k_j$ for all $i, j$ and denoting $T_{12}:=T_{1down} =T_{2up}$, we have $$dot{T_{1}} = {k} (T_1 - frac{T_{1up}+ T_{1down}}{2}) approx {k} (T_1 - frac{T_{1up}+ T_{12}}{2}),$$ and similarly $$dot{T_{2}} approx {k} (T_{2} - frac{T_{12}+ T_{4}}{2}),$$ $$dot{T_{3}} approx {k} (T_3 - frac{T_{4}+ T_{3down}}{2}),$$ so that the approximation $T_4 approx frac{T_2 + T_3}{2}$ implies that $$dot{T_4} = k (frac{T_4}{2} - frac{T_{12} + 2 T_4 + T_{3down}}{4}) = - k frac{T_{1} + T_{2} + 2 T_{3down}}{8} < 0.$$ This implies that $dot{T_4} rightarrow 0$ as $t rightarrow 0$ if $T_2$ and $T_3$ are stable. Considering that temperature PID control of zones $2$ and $3$ is available, a good first approach is to apply individual PID controllers for zones $1$ and $2$ to stabilize their temperatures and set $T_2^{ref}$ and $T_3^{ref}$ such that $T_4^{ref} = frac{T_2^{ref} + T_3^{ref}}{2}$.