I'm trying to recreate Figure 1 in this paper. This requires maximizing equation (19), which I have convinced myself is concave, but I am having trouble implementing it in CVXPY. Here is the code I think is relevant to this question:
import cvxpy as cp import numpy as np def tanh_atom(x): return (cp.exp(2 * x) - 1) / (cp.exp(2 * x) + 1) def functional(X, T, dt, G=None): if G is None: G = np.identity(X.shape-1) Q = np.linalg.inv(np.matmul(G, G.T)) dXk_dTk = cp.diff(X, axis=0) / dt Fk = tanh_atom(X)[:-1] diff = dXk_dTk - Fk S = dt * cp.quad_form(diff, Q) return -.5 * S def norm_pdf_atom(x, mu, sd): coef = 1 / (sd * (2*np.pi)**.5) z = (x - mu) / sd return coef * cp.exp(-.5 * cp.power(z, 2)) def merit(obs, X, T, dt): term1 = functional(X, T, dt) term2 = norm_pdf_atom(X, 0, .16**.5) term3 = norm_pdf_atom(obs, X[-1], .16**.5) return term1+term2+term3 t0 = 0 t_end = 5 N = 5 obs = 1.5 T = np.reshape(np.linspace(t0, t_end, N), (N, 1)) dt = T[:,0] - T[:,0] dt = cp.Parameter(value=dt, nonneg=True) T = cp.Parameter(np.shape(T), value=T, nonneg=True) X = cp.Variable(np.shape(T)) constraints = [X >= 0, X <= 10] prob = cp.Problem(cp.Maximize(merit(obs,X,T,dt)), constraints) prob.solve(verbose=True)
- CVXPY tells me that
term1is not DCP siting this expression:
(exp(Promote(2.0, (5, 1)) @ var2) + Promote(-1.0, (5, 1))) / (exp(Promote(2.0, (5, 1)) @ var2) + Promote(1.0, (5, 1)))
- For both
term3CVXPY tells me that the expression below is not DCP. It actually suggests using DQCP with
solve()call but this also gives an error. Both these terms are just maximizing a Gaussian distribution so I really feel like this shouldn't be an issue. Is there a better way to do this other than creating the pdf function?
exp(-0.5 @ power((var2[0, 0] + -0.0) / 0.4, 2.0))
In both scenarios I can't figure out why they are not DCP, as they both are concave, unless I messed up my reasoning for that. I've tried using as many CVXPY atoms as possible but still don't know what's going on. Any help would be much appreciated.
Note: I individually tested each term in
merit just by commenting out the other terms and running the script to see if there were errors. I have also recreated this situation using
numpy by changing every place in the script with
np. Then I would randomly generate column vectors
X and feed them to
merit and plot the best
X. Although this is a very crude way to find the optimal
X, it gave me results close enough to the desired figure to make me think the math behind the code is correct.
Edit: I found some errors in my original question so I have made changes as necessary.