You can create a complex, nonlinear, numerical model of the full system in as much detail as you want (eg implemented in a block sim like Simulink, or MBD software like ADAMS, or just A and B matrices derived manually but incorporating nonlinearities like sin, etc).
From that you can evaluate a whole bunch of finite differences around an arbitrary state (linearization point), which are the components of the A and B matrices. This is an approximate, linear model which is "close enough" in behaviour to your real model for states that are "close" to your linearization point.
In practice if you're moving far away from your linearization point (eg a 6 DOF robot arm with large displacement of all joints), you'll need to re-compute the linearization each timestep, around the current state, and re-compute the LQR or MPC solution from this "new" approximate model.
Here's some linearization code I wrote for a library that's used to teach controls and robotics in undergrad. It takes as a input a "model" that's an arbitrary function that generates the next state from a given state, and computes the jacobian matrices by finite differences: https://github.com/SherbyRobotics/pyro/blob/43dcb112427978ff...