Issue solving nonlinear coupled PDEs and optimization

I am trying to solve a set of coupled PDEs (equations in attached image), where the dependent variables need to satisfy a few constraints, and the boundary conditions change with the solution of PDEs with time (optimization problem). I think I need to solve incrementally (in time) using pdpe and check for constraints, update my boundary conditions, and pass the last array values as new initial conditions for next step. However, I do not have an equation for explicitly, so will that be a problem if i am to use pdpe solver? If not, what other methods are feasible?

Edit: Can I discretise the equations so i have non linear set of algebraic equations, which I can solve using fsolve? I have never seen fsolve be used for such a case, so would it be feasible to try it? I would have around 300 variables, for 62 grid points.

[Torsten]

What is the equation for "phi" ?

[Saketharam Reddy N]

I do not have an explicit equation for ''phi''. The relation is implicit as given in the equations i posted.

Edit: I may be using the terms implicit and explicit incorrectly, but what i want to say is that there are no other equations apart from the ones i have posted.