Qutip 4 feedback solver

qutip4_feedback.py

## From https://github.com/qutip/qutip/issues/1571#issuecomment-859873615

times = [0, ...]
# Liouvillians for parts you have control over
control_liouvillians = [qutip.liouvillian(h) for h in control_hamiltonians]
# The time-dependent Hamiltonian for stuff you're not controlling
base_hamiltonian = qutip.QobjEvo([H0, [H1, time_dependence], ...])
# Turn it into a Liouvillian once, so we don't repeat the cost
base = qutip.liouvillian(base_hamiltonian, collapse_operators)
state = ...
options = qutip.Options(store_states=False, store_final_state=True)
for prev, time in zip(times[:-1], times[1:]):
    controls = krotov.get_next_controls(time, state, ...)
    current_liouvillian = base.copy()
    for control, operator in zip(controls, control_liouvillians):
        current_liouvillian += control * operator
        #                      ^^^^^^^^^^^^^^^^^^
        # each of these terms is a single time-independent Qobj,
        # and the sum is a single QobjEvo with all the uncontrolled
        # time dependence already handled.
    state = qutip.mesolve(current_liouvillian, state, [prev, time], options=options).final_state