Maintaining the stability and shape of a plasma is a crucial task in many technological applications ranging from beam shaping to fusion energy. This is often challenging as plasma systems tend to be naturally unstable and kinetic effects can play an important role in the behavior of the instabilities. Due to the large energies involved, the primary way to control plasma systems is the application of external electromagnetic fields.In this work, we deploy a PDE-constrained optimization formulation that uses a kinetic description for plasma dynamics as the constraint. This is to optimize, over all possible controllable external electric fields, the stability of the plasma dynamics under the condition that the Vlasov–Poisson (VP) equation is satisfied. Mathematically we formulate it as a PDE-constrained optimization. For computing the functional derivative with respect to the external field in the optimization updates, the adjoint equation is derived. Furthermore, in the discrete setting, where we employ the semi-Lagrangian method as the forward solver, we also explicitly formulate the corresponding adjoint solver and the gradient as the discrete analogy to the adjoint equation and the Fréchet derivative. A distinct feature we observed of this constrained optimization is the complex landscape of the objective function and the existence of numerous local minima, largely due to the hyperbolic nature of the VP system. To overcome this issue, we utilize a gradient-accelerated genetic algorithm, leveraging the advantages of the genetic algorithm's exploration feature to cover a broader search of the solution space and the fast local convergence aided by the gradient information. We show that our algorithm obtains good electric fields that are able to maintain a prescribed profile in a beam shaping problem and uses nonlinear effects to suppress plasma instability in a two-stream configuration.