  We present a finite-horizon optimization algorithm that extends the
established concept of Dual Dynamic Programming (DDP) in two ways. First, in
contrast to the linear costs, dynamics, and constraints of standard DDP, we
consider problems in which all of these can be polynomial functions. Second, we
allow the state trajectory to be described by probability distributions rather
than point values, and return approximate value functions fitted to these. The
algorithm is in part an adaptation of sum-of-squares techniques used in the
approximate dynamic programming literature. It alternates between a forward
simulation through the horizon, in which the moments of the state distribution
are propagated through a succession of single-stage problems, and a backward
recursion, in which a new polynomial function is derived for each stage using
the moments of the state as fixed data. The value function approximation
returned for a given stage is the point-wise maximum of all polynomials derived
for that stage. This contrasts with the piecewise affine functions derived in
conventional DDP. We prove key convergence properties of the new algorithm, and
validate it in simulation on two case studies related to the optimal operation
of energy storage devices with nonlinear characteristics. The first is a small
borehole storage problem, for which multiple value function approximations can
be compared. The second is a larger problem, for which conventional discretized
dynamic programming is intractable.
