A method of controlling a system that has a finite execution time in the optimization of a cost function: J = ∑ i = 1 N ⁢ F ⁡ ( u ⁡ ( k + i | k ) , y ⁡ ( k + i | k ) ) , where N is selected such that J is convergent; k+i is the time in the future after the passage of i time increments; u(k+i|k) and y(k+i|k) are restricted to values of the respective control variables u and performance variables y that are compatible with predetermined system operational constraints; and for each time increment i=1 . . . N, values for F(u(k+i|k), y(k+i|k)) are calculated using a system model M to determine y(k+i|k) for a plurality of respective discrete values of u(k+i|k) that are limited to and representatively sampling the range of possible values of u that are compatible with the predetermined system operational constraints, J being optimised by identifying the combination {u(k+1|k),u(k+2|k), . . . u(k+N|k)} of respective discrete values u(k+i|k) over N time increments that combination yields an optimised value for J.