From d216b754b2e66fe33abfa83671cb0446c225b104 Mon Sep 17 00:00:00 2001 From: Ethan Lew Date: Sat, 13 Apr 2024 12:46:17 -0700 Subject: [PATCH] add regularization because of sparse case (@Abdu-Hekal #88) --- autokoopman/estimator/koopman.py | 15 ++++++++++++--- notebooks/weighted-cost-func.ipynb | 14 +++++++------- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/autokoopman/estimator/koopman.py b/autokoopman/estimator/koopman.py index e8dbeb3..4507f7b 100644 --- a/autokoopman/estimator/koopman.py +++ b/autokoopman/estimator/koopman.py @@ -55,6 +55,8 @@ def wdmdc(X, Xp, U, r, W): def swdmdc(X, Xp, U, r, Js, W): """State Weighted Dynamic Mode Decomposition with Control (wDMDC)""" import cvxpy as cp + import warnings + assert len(W.shape) == 2, "weights must be 2D for snapshot x state" if U is not None: @@ -76,7 +78,8 @@ def swdmdc(X, Xp, U, r, Js, W): # SW-eDMD objective weights_obj = np.vstack([(np.abs(J) @ w)[:n_obs] for J, w in zip(Js, W)]).T P = sf * cp.multiply(weights_obj, Xp[:, :n_obs].T - K @ X[:, :n_obs].T) - objective = cp.Minimize(cp.sum_squares(P)) + # add regularization + objective = cp.Minimize(cp.sum_squares(P) + 1E-4 * 1.0 / (n_obs**2) * cp.norm(K, "fro")) # unconstrained problem constraints = None @@ -85,8 +88,14 @@ def swdmdc(X, Xp, U, r, Js, W): prob = cp.Problem(objective, constraints) # solve for the SW-eDMD Koopman operator - result = prob.solve() - # TODO: check the result + _ = prob.solve(solver=cp.CLARABEL) + #_ = prob.solve(solver=cp.ECOS) + + # backup case + if K.value is None: + # give a warning about the optimization failure + warnings.warn("SW-eDMD (cvxpy) Optimization failed to converge. Switching to unweighted DMDc.") + return dmdc(X, Xp, U, r) # get the transformation Atilde = K.value diff --git a/notebooks/weighted-cost-func.ipynb b/notebooks/weighted-cost-func.ipynb index 8ea8c0a..6a22183 100644 --- a/notebooks/weighted-cost-func.ipynb +++ b/notebooks/weighted-cost-func.ipynb @@ -45,7 +45,7 @@ "import autokoopman.benchmark.fhn as fhn\n", "fhn = fhn.FitzHughNagumo()\n", "training_data = fhn.solve_ivps(\n", - " initial_states=np.random.uniform(low=-2.0, high=2.0, size=(10, 2)),\n", + " initial_states=np.random.uniform(low=-2.0, high=2.0, size=(1, 2)),\n", " tspan=[0.0, 6.0],\n", " sampling_period=0.1\n", ")\n", @@ -90,7 +90,7 @@ " #w = np.sum(traj.abs().states, axis=1)\n", " #w = 1/(traj.abs().states+1.0)\n", " w = np.ones(traj.states.shape)\n", - " w[:, -3:] = 1.0\n", + " w[:, -3:] = 0.0\n", " w[:, :2] = 1.0\n", " weights.append(w)\n", "\n", @@ -136,8 +136,8 @@ " n_obs=40, # maximum number of observables to try\n", " max_opt_iter=200, # maximum number of optimization iterations\n", " grid_param_slices=5, # for grid search, number of slices for each parameter\n", - " n_splits=None, # k-folds validation for tuning, helps stabilize the scoring\n", - " rank=(20, 21, 1) # rank range (start, stop, step) DMD hyperparameter\n", + " n_splits=None, # k-folds validation for tuning, helps stabilize the scoring\n", + " rank=(40, 41, 1) # NOTE: don't rank tune (for now)--SW-eDMD already optimizes for this\n", ")" ] }, @@ -168,9 +168,9 @@ " n_obs=40, # maximum number of observables to try\n", " max_opt_iter=200, # maximum number of optimization iterations\n", " grid_param_slices=5, # for grid search, number of slices for each parameter\n", - " n_splits=5, # k-folds validation for tuning, helps stabilize the scoring\n", + " n_splits=None, # k-folds validation for tuning, helps stabilize the scoring\n", " lengthscale=(0.1, 1.0),\n", - " rank=(20, 21, 1) # rank range (start, stop, step) DMD hyperparameter\n", + " rank=(1, 41, 5) # rank range (start, stop, step) DMD hyperparameter\n", ")" ] }, @@ -290,7 +290,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.9.0" } }, "nbformat": 4,