I solve the problem of optimal control. First, I solve a system of two equations (with respect to the variables p, q, V, n) with already optimal conditions (Kopt and Sopt) to check whether it is possible to solve such a system at all. It turns out that everything is OK, but the non-linear part of the equation cannot be discarded, it is important. I am using the NDSolve function with the StiffnessSwitching specification:

`eq1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0; eq2 = D[V[x, t], t] + Kopt[t]*x - 2*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] + g*D[n[x, t], x] - D[ D[Sqrt[n[x, t]], x, x]/(2*Sqrt[n[x, t]]), x] == 0; eqo1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0; eqo2nl = D[V[x, t], t] + Kopt[t]*x - 8*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] + g*D[n[x, t], x] - ((D[n[x, t], x])^3 - 2*n[x, t]*D[n[x, t], x]*D[n[x, t], {x, 2}] + D[n[x, t], {x, 3}]*(n[x, t])^2)/(4*(n[x, t])^3) == 0; cond1 = n[x, t] == ((Sqrt[Kzero]/Pi)^(1/2))* Exp[-0.5*Sqrt[Kzero]*x^2] /. t -> 0; cond2 = n[x, t] == 0. /. x -> -10; cond3 = n[x, t] == 0. /. x -> 10; cond4 = V[x, t] == 0. /. t -> 0; cond5 = V[x, t] == 0. /. x -> -10; cond6 = V[x, t] == 0. /. x -> 10; sol5 = NDSolve[{eqo1, eqo2nl, cond1, cond2, cond3, cond4, cond5, cond6}, {V, n}, {x, -10, 10}, {t, 0, Tv}, Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 1, PrecisionGoal -> 1]; `

Next, I solve four equations in order to find the control parameter K (depends on q), while Sopt is still known. The initial conditions are the same and as banal as possible. Add two more equations and K[t] formula:

`K[t_] := -1*Integrate[q[x, t]*x^2, {x, -10, 10}]; eqo3 = D[q[x, t], t] == - n[x, t]*D[p[x, t], x] - V[x, t]*D[q[x, t], x]; eqo4 = D[p[x, t], t] + D[p[x, t], x]*V[x, t] + g*D[q[x, t], x] == 0; `

The linearized system has a solution, but not a nonlinear one! There are many problems, including the order of the equation:

`The spatial derivative order of the PDE may not exceed two. `

I try to remove the third order derivative from the nonlinear part:

`The maximum derivative order of the nonlinear PDE coefficients for \ the Finite Element Method is larger than 1. It may help to rewrite \ the PDE in inactive form. `

Why does he even use only FEM now? Why was a nonlinear system of two equations solved normally? What should I do, help me, please.