[Baseline] Linear optimal control for cylinder

[jcallaham]

Goal is to have a working LQG (LQR + Kalman filter) for the cylinder wake as a baseline controller. So far I've gotten the following to work (see examples/cylinder/control and examples/cylinder/notebooks/controller-design.ipynb):

• Convert linearized timestepper to a discrete-time LTI system with the matrix-free LinearOperator from scipy
• Extract direct and adjoint global eigenmodes from SLEPc as numpy arrays
• Derive a reduced-order model with Petrov-Galerkin projection onto the global modes
• Controller design with the Python Control Systems Library
• Tested both the Kalman filter works and full LQG controller on the ROM

So far so good, except that it blows up when the controller is actually applied to the DNS. Here's my plan to isolate the breakdown:

• [x] Test Kalman filter and LQR on the ROM as a sanity check
• [x] Redo everything with the smaller mesh so things can be tested more efficiently
• [x] Kalman filter on full LTI system
• [ ] LQR on full LTI system
• [ ] Kalman filter on full linearized timestepper (not LinearOperator LTI)
• [ ] LQR on full linearized timestepper
• [ ] Kalman filter on nonlinear timestepper
• [ ] LQR on nonlinear timestepper

[jcallaham]

Started redoing with the smaller mesh and ran into a new issue that seems to have to do with parallelization in extracting the eigenvectors from SLEPc. Made kind of a mess of it so far, but I think I found something that might be helpful for Modred as well (#13). Basically I think we can "scatter" the PETSc.Vec object to rank zero, and then save as a numpy binary.

• https://petsc.org/main/docs/manualpages/PetscSF/VecScatterCreateToZero/
• https://github.com/erdc/petsc4py/blob/e29612c871aa29c8193c7559bbe425f4ce920bbc/src/PETSc/Scatter.pyx#L51-L57

[jcallaham]

Turns out this is a little more complicated because of the local/global mappings... haven't been able to figure out how to scatter to rank zero with the same ordering that a serial run would have. But maybe this isn't that important here... there are really two issues.

1. Extracting lower-level data to binary for switching from complex to real mode
2. Gathering data from parallel SLEPc analysis to a firedrake Function

I've basically been trying to do them both in one shot, but maybe the easiest thing to do for now would be to do the real/imag split with the vectors split into parallel (work with q.dat.data or vec.array, for instance). Then save from each rank, reload in "real mode", rebuild the PETSc vector, and assign to the firedrake Function.

[jcallaham]

This seems to work, at least for loading to/from a Function in parallel. Since the intermediate is a numpy array, no reason why this also wouldn't work to go from complex to real.

q1 = flow.q.copy(deepcopy=True)
norm1 = flow.dot(ufl.real(q1), ufl.real(q1))

with q1.dat.vec_ro as vec:
    rank = fd.COMM_WORLD.rank
    np.save(f'tmp/vec_{rank}.npy', vec.array)

# And reload
q2 = fd.Function(flow.mixed_space)
with q2.dat.vec as vec:
    rank = fd.COMM_WORLD.rank
    vec.array = np.load(f'tmp/vec_{rank}.npy')

norm2 = flow.dot(ufl.real(q2), ufl.real(q2))
gym.print(np.isclose(norm1, norm2))

I think the thing to do would be to add a utility function like gather that can read in the rank-split array and rebuild a Function. Then maybe stability analysis would just run with a bash script so that the same number of processors is guaranteed to be used for both.

[jcallaham]

This worked! And I think is a much better solution than the original anyway. Still need to make a utility function to clean that up, but it's good for the time being. Moving on to testing the full-order LTI system with estimation/control

[jcallaham]

Confirmed the Kalman filter works with linearized LTI system, but adding control causes a blow-up:

Not sure if the issue is the controller itself or the implementation in the simulation. One thing to try might be using the step functionality rather than the A*x + B*u form with LinearOperators. Or some kind of limiting or smoothing in the time-stepping actuation model.

[jcallaham]

Also check if A*x + B*u form is equivalent to the step functionality with linearized time-stepper...

[jcallaham]

Now I think the issue might be that the control is so tightly coupled to the measurement here, since the actuation and measurement both just act on the surface of the cylinder... not sure how to handle that exactly. This might even just be a numerical issue since the timestepper has to respond to the strongly enforced Dirichlet BCs.

One possibility would be to return observations as a moving average, similar to how control is done.