  We propose an algorithm to solve optimization problems constrained by partial
(ordinary) differential equations under uncertainty, with almost sure
constraints on the state variable. To alleviate the computational burden of
high-dimensional random variables, we approximate all random fields by the
tensor-train decomposition. To enable efficient tensor-train approximation of
the state constraints, the latter are handled using the Moreau-Yosida penalty,
with an additional smoothing of the positive part (plus/ReLU) function by a
softplus function. In a special case of a quadratic cost minimization
constrained by linear elliptic partial differential equations, and some
additional constraint qualification, we prove strong convergence of the
regularized solution to the optimal control. This result also proposes a
practical recipe for selecting the smoothing parameter as a function of the
penalty parameter. We develop a second order Newton type method with a fast
matrix-free action of the approximate Hessian to solve the smoothed
Moreau-Yosida problem. This algorithm is tested on benchmark elliptic problems
with random coefficients, optimization problems constrained by random elliptic
variational inequalities, and a real-world epidemiological model with 20 random
variables. These examples demonstrate mild (at most polynomial) scaling with
respect to the dimension and regularization parameters.
