ipopt report

7/report7.md

@@ -23,24 +23,26 @@ let `u` = control state on boundary, excl. corners
          1. We get one equation per **inner** grid point, $G_{ij}=0$
       1. The control $y=u~\text{on}~\partial\Omega$
          1. We get this $4n$ times 
-         1. $\implies$ total number of constraints $M = n^2 + 4n < N = |x|$, so there are **less constraints than variables!**
+         1. $\implies$ total number of constraints $M = n^2 + 4n < N = |x|$, so there are **less constraints than variables!** This means there are multiple feasible solutions... of which we want to find the optimal one.
    1. The `jacobian` (of the constraints) is pretty complicated. Jacobian sparsity pattern <img src="jacobian-sparsity.png" style="zoom: 30%;" />
    1. The `hessian` (of the lagrangian) is just a $N \times N$ diagonal
-1. Choose initial state vector to be *anything that satisfies all the constraints*. Starting with $y^d$ this became $y_0(x_1, x_2) = -10x_1(x_1 - 1)x_"(x_2 - 1) + 1$. With the initial $u$ being equal to the $y_0$ on the boundary (to satisfy the constraints).
-   1. Not sure about the sign! The initial conditions have a large impact on the solution. Different initial states generally do not converge to the same solution.
-   1. Positive sign (top=$y_0$, bottom=converged solution):  
+1. Choose initial state vector to be anything where the result is feasible. A state that satisfies all the constraints makes sense, but is not necessary. It will however likely require less iterations to converge, so is preferred.  
+   For me, starting with $y^d$ this became $y_0(x_1, x_2) = -5x_1(x_1 - 1)x_"(x_2 - 1) + 1$. With the initial $u$ being equal to the $y_0$ on the boundary (to satisfy the control constraints).
+   1. You should check whether the constraints are satisfied at the end. In this case that they should all evaluate to roughly 0. In the code this is checked by calculating the norm of the constraint vector `constraint_norm`.
+   1. The initial conditions have a large impact on the solution. Different initial states generally do not converge to the same solution.
+   1. (top=$y_0$, bottom=converged solution):  
       ![](positive-y0.png)
-1. Negative sign:  
-       ![](negative-y0.PNG)   
-      (We can see that the constraint $y(x) < 3.5$ is working)
+   1. Different $y_0$ with negative sign:  
+          ![](negative-y0.PNG)   
+         (We can see that the constraint $y(x) < 3.5$ is working)
 1. See plots, the purple points on the edges denote $u$, the surface denotes $y(x)$
 1. Both the hessian of the lagrangian and the jacobian of the constraints are constant. We can notify IPOPT by using the `nlp.addOption('hessian_constant', 'yes')` and `nlp.addOption('jac_c_constant', 'yes')`. These will then only be evaluated once.
-1. Making both the hessian and jacobian *sparse* made the largest difference in performance, as these matrices are very large and sparse. This is done in a triplet format using `hessianstructure()` and `jacobianstructure()`.
+1. Making both the hessian and jacobian *sparse* made the largest difference in performance, as these matrices are very large and sparse. This is done in a triplet format using `hessianstructure()` and `jacobianstructure()`. It is also important that the construction does not rely on dense matrices and `np.nonzero`.
    1. Note that the size of the matrices are in $\mathcal{O}(n^4)$
-   1. It seems as is there is no effect of adding more cores. However, the variance is very low so there is some difference. Pretty sure there is something wrong here.![](perf-results.png)
-   1. *Most of the time (>50%) is spent initializing*, the time for the iterations is relatively low
-   1. Testing grid sizes of n=500-3000 is infeasible with this implementation
-1. Short answer: Its easy to enforce the pde constraint when including $y$.
+   1. It seems as is there is little effect of adding more cores. The variance is generally very low, so only one repeat is shown here.![](perf-results.png)
+   1. Initializing also takes some time, e.g. for $n=1000$ we have ~4s.
+1. Short answer: Its easy to enforce the pde constraint when including $y$. Faster to solve? Other library is needed otherwise.
+   1. Alternative: reduce-space instead of full-space approach, requires a quasi-Newton method
    1. This would make the jacobian and hessian much smaller, also asymptotically they will be $\mathcal{O}(n^2)$ as the size of the state `x` is then $\mathcal{O}(n)$.  
       The hessian would still be a diagonal with $h\alpha$.  
       The jacobian would be a diagonal of ones? These would give large performance benefits.