1
$\begingroup$

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[0]-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][0], 0, .16**.5)
    term3 = norm_pdf_atom(obs, X[-1][0], .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][1] - T[:,0][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)

Two issues

  • CVXPY tells me that term1 is 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 term2 and term3 CVXPY tells me that the expression below is not DCP. It actually suggests using DQCP with qcp=True in the 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) / 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 cp to 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.

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy