A stochastic model predictive controller for controlling a system under uncertainty subject to constraints on state and control variables of the system, comprising: at least one processor; and memory having instructions stored thereon that, when executed by the at least one processor, cause the stochastic model predictive controller to: solve, at each control step of the stochastic model predictive controller, an inequality constrained nonlinear dynamic optimization problem including probabilistic chance constraints representing the uncertainty to produce a control command, wherein the constraints on state and control variables of the system include inequality constraints, wherein the probabilistic chance constraints are based on a tightening of one or multiple inequality constraints with a term that depends on a back-off coefficient value ( α i 420), a constraint Jacobian matrix and a covariance matrix for the predicted state variables at each time step within the prediction horizon, to ensure that the probability of violating each of the corresponding inequality constraints is below a probability threshold value ( ε i), and wherein the dynamic optimization problem is solved based on a two-level optimization that alternates, until a termination condition is met, propagation of the covariance matrices of the probabilistic chance constraints at each time step within the prediction horizon for fixed values of the state and control variables with optimization of the state and control variables within the prediction horizon for fixed values of the covariance matrices; and control an operation of the system using the produced control command including the optimized control variables, wherein the propagation of the covariance matrices for each time step within the prediction horizon is performed by evaluating nonlinear covariance propagation equations with respect to current values of the state and control variables optimized for the prediction horizon, wherein the nonlinear covariance propagation equations are linearization-based discrete-time Lyapunov equations that preserve positive definiteness of the covariance matrices at each control time step, using Jacobian matrices precomputed with respect to the current values of the state and control variables optimized for the prediction horizon. solve, at each control step of the stochastic model predictive controller, an inequality constrained nonlinear dynamic optimization problem including probabilistic chance constraints representing the uncertainty to produce a control command, wherein the constraints on state and control variables of the system include inequality constraints, wherein the probabilistic chance constraints are based on a tightening of one or multiple inequality constraints with a term that depends on a back-off coefficient value ( α i 420), a constraint Jacobian matrix and a covariance matrix for the predicted state variables at each time step within the prediction horizon, to ensure that the probability of violating each of the corresponding inequality constraints is below a probability threshold value ( ε i), and wherein the dynamic optimization problem is solved based on a two-level optimization that alternates, until a termination condition is met, propagation of the covariance matrices of the probabilistic chance constraints at each time step within the prediction horizon for fixed values of the state and control variables with optimization of the state and control variables within the prediction horizon for fixed values of the covariance matrices; and control an operation of the system using the produced control command including the optimized control variables, wherein the propagation of the covariance matrices for each time step within the prediction horizon is performed by evaluating nonlinear covariance propagation equations with respect to current values of the state and control variables optimized for the prediction horizon, wherein the nonlinear covariance propagation equations are linearization-based discrete-time Lyapunov equations that preserve positive definiteness of the covariance matrices at each control time step, using Jacobian matrices precomputed with respect to the current values of the state and control variables optimized for the prediction horizon.