The joint moments of the derivatives of the characteristic polynomial of a random unitary matrix, and also a variant of the characteristic polynomial that is real on the unit circle, in the large matrix size limit, have been studied intensively in the past twenty five years, partly in relation to conjectural connections to the Riemann zeta-function and Hardy's function. We completely settle the most general version of the problem of convergence of these joint moments, after they are suitably rescaled, for an arbitrary number of derivatives and with arbitrary positive real exponents. Our approach relies on a hidden, higher-order exchangeable structure, that of an exchangeable array. Using these probabilistic techniques, we then give a combinatorial formula for the leading order coefficient in the asymptotics of the joint moments, when the power on the characteristic polynomial itself is a positive real number and the exponents of the derivatives are integers, in terms of a finite number of finite-dimensional integrals which are explicitly computable. Finally, we develop a method, based on a class of Hankel determinants shifted by partitions, that allows us to give an exact representation of all these joint moments, for finite matrix size, in terms of derivatives of Painlevé V transcendents, and then for the leading order coefficient in the large-matrix limit in terms of derivatives of solutions of the $σ$-Painlevé III' equation. Equivalently, we can represent all the joint moments of power sum linear statistics of a certain determinantal point process behind this problem in terms of derivatives of $σ$-Painlevé III' transcendents. This gives an efficient way to compute all these quantities explicitly. Our methods can be used to obtain analogous results for a number of other models sharing the same features.