Rigid robots powered by electric or hydraulic actuation are often bulky and unsuitable for propulsion in soft, uneven substrates. In contrast, magnetically actuated soft robots offer a new class of tether-free locomotion in complex terrains and fluid environments. These systems rely on magnetoelastic materials such as soft silica gels embedded with magnetizable oxides that exhibit highly nonlinear, multiphysics dynamics governed by elasticity and electromagnetism [1]. Incorporating such dynamics into a control framework naturally leads to partial differential equation (PDE) constrained optimization challenges. Despite significant advances in the modelling of soft robotic dynamics, several important gaps remain unexplored in their control. In particular, the conservation laws associated with the symmetries of magnetoelastic control systems are still not well understood. Likewise, numerical discretization for the optimal control of magnetoelastic PDEs that preserves the geometric invariants of these systems is not yet available [2, 3]. From a computational standpoint, control of such systems remains challenging due to their multiphysics coupling and the high dimensionality of the resulting nonlinear programming problems (NLPs). Thus, the integration of PDE-constrained optimization with structure-preserving algorithms, ensuring that conservation laws are respected during control, has not been systematically addressed. In this study, we present an adjoint-based optimal control framework for magnetoelastic thin, slender structures. The external magnetic field is treated as the control variable, enabling contactless actuation of the system. Conservation laws associated with the symmetries of magnetoelastic control systems are deduced through Noether’s theorem [2, 4]. The governing PDEs are discretized in space using finite elements, resulting in time-continuous differential-algebraic equations (DAEs). Next, we derive the necessary optimality conditions for a tracking-type objective functional. Finally, to ensure numerical stability, we employ structure-preserving integration of the optimality conditions and solve the resulting NLP efficiently using the forward-backward sweep method (FBSM) [4, 5].