You are dealing with a multi-input multi-output system (MIMO) system in which case you can't in general apply the same techniques as for single-input single-output (SISO) systems. If you are dealing with a linear time invariant MIMO system, then you could use something called the relative gain array to see how well a SISO control approach would work. However in this case I believe you are dealing with a nonlinear system, but it might be possible to do input-output decoupling, namely by controlling the netto inflow into each tank. For this I assume for now that you can choose the flow rates $q_0$, $q_1$, $q_2$ and $q_3$ and that the dynamics can be written as
$$
\begin{align}
\dot{h}_i = \alpha_i (q_{i-1} - q_i) \tag{1}
\end{align}
$$
so the combines dynamics can be written as
$$
\dot{\vec{h}} =
\begin{bmatrix}
\alpha_1 & -\alpha_1 & 0 & 0 \\
0 & \alpha_2 & -\alpha_2 & 0 \\
0 & 0 & \alpha_3 & -\alpha_3
\end{bmatrix} \vec{q} \tag{2}
$$
where $\vec{h} = \begin{bmatrix}h_1 & h_2 & h_3\end{bmatrix}^\top$ and $\vec{q} = \begin{bmatrix}q_0 & q_1 & q_2 & q_3\end{bmatrix}^\top$. Now by defining a virtual input $\vec{v}$ such that $\dot{\vec{h}} = \vec{v}$, then every input and output pair is decoupled. For this we need
$$
\begin{bmatrix}
\alpha_1 & -\alpha_1 & 0 & 0 \\
0 & \alpha_2 & -\alpha_2 & 0 \\
0 & 0 & \alpha_3 & -\alpha_3
\end{bmatrix} \vec{q} = \vec{v} \tag{3}
$$
which can be solved for $\vec{q}$ using a pseudo inverse of the matrix on the left hand side. There are multiple ways this could be done, but the solution which minimizes $\|\vec{q}\|$ yields
$$
\vec{q} =
\begin{bmatrix}
\frac{3}{4\,\alpha_1} & \frac{1}{2\,\alpha_2} & \frac{1}{4\,\alpha_3} \\
\frac{-1}{4\,\alpha_1} & \frac{1}{2\,\alpha_2} & \frac{1}{4\,\alpha_3} \\
\frac{-1}{4\,\alpha_1} & \frac{-1}{2\,\alpha_2} & \frac{1}{4\,\alpha_3} \\
\frac{-1}{4\,\alpha_1} & \frac{-1}{2\,\alpha_2} & \frac{-3}{4\,\alpha_3}
\end{bmatrix} \vec{v}. \tag{4}
$$
So now you can use normal SISO control methods, like PID, for each $h_i$-$v_i$ pair which would drive $\vec{h}$ to a desired value.
However there is probably not a linear relation between $q_i$ and $V_i$, but more likely something like $q_i \propto \sqrt{V_i(h_i - h_{i+1})}$ or some other nonlinear state depended function. If you know this relation more accurately you could come up with a nonlinear controller, otherwise you could perform sequential loop closing. For the sequential loop closing you use $V_i$ to control $q_i$ to its desired value defined by $\vec{v}$ from equation $(4)$.
If you have some more challenging reference signal for the heights, then you might have to resort to even more advance techniques, such as model predictive control. For example when the reference heights are such that $h_{r,1} < h_{r,2}$ and $h_{r,2} > h_{r,3}$ but your system starts at $h_1 = h_2 = h_3 < h_{r,2}$, then it will initially not be possible to have a positive inflow into the second tank and to reach the desired heights $h_1$ actually has to be bigger than $h_2$ for a while in order for $h_2$ to reach its goal.
