Some responses:

"Your code is not using dsolve numeric in parametric form. I will be posting the use of parametric form working with NLPSolve (thanks to Allan Wittkopf)."

I agree - as I said, I thought from your original post you may have been looking for a workaround.

As I said, we don't want analytical solutions.

This was just to test my answer.

Also, one parameter optimization is different from multiple parameters (scalar vs vector).

Sorry, I have misunderstood the problem statement. In your reply to @vv you said 

"The problem statement is given below

dca/dt = -(u+u^2/2)*ca

dcb/dt = u*ca - 0.1*cb

t = 0, ca = 1, cb = 0
Maximize cb(1) subject to the bounds 0 < u< 5."

which seemed to me to have only one parameter.

My method is easily extended to more parameters.

To be clear, I am not looking for an answer for the model. My goal was to use NLPSolve with dsolve/numeric/parameters with sqp.

Yes, I understood that was just a test problem, and I was proposing a generic solution. In view of your comments about the complexity of the solution for least squares, I stripped done my more general earlier answer to make it more clear how it works.

But anyway it seems you now have a solution. Good luck!

@ The code is that way because I had it already at hand, and I didn't have a lot of time this morning. The complexity of the code was to encapsulate the whole thing in a user-friendly callable routine without using global variables. I realize that was not your issue, but since you were having problems with the dsolve parametric form, I thought you might be interested in this alternative way of doing it, which just substitutes the parameters in.

It is still easier for me to rewrite than debug your code. I'm guessing somehow the nvars=2 problem just doesn't have enough degrees of freedom to be sensibly done by NLPSolve, since it works with nvars higher. Anyway, here is a solution to your problem if I have understood it correctly. If you need the function cb at various specified times, you can use the output=Array option in the final dsolve call.

I went back to the procedure form, since then you can let the dsolve error mechanism decide on the stepsize etc and bypass using nvar. In my experience the Matrix method isn't much faster, but that depends on the class of problems.

Download defit5.mw

I now tried it with sqp, and it seems to work, but is slower.

@MaPal93 It turns out that it is a solution, one just needs more accuracy to show that:

Numerical_issue.mw

The small change in length I think just means Maple generated and stored the expression in a slightly different form.

@MaPal93 The rho_u case seems to have a positive solution, but when I substitute into EqN_rhou I don't find it is a solution. Perhaps the solution matches another version of EqN_rhou?

250523_SOLUTION_analytical2.mw

@MaPal93

"Does it mean that the quadratic solutions I obtain from the solver would be the same regardless of the assumptions on my variables and parameters (and only the roots to these quadratic equations depend on the assumptions)?" 

Yes, I do not think PolynomialSystem is using the assumptions and so you could (should) leave them out. In principle if you put them in, say by using RegularChains:-Triangularize directly (where you can just add them to the set of equations), and it could figure out that the quadratic did not have any positive solutions then it would not return the quadratic as a solution. I only tried RegularChains:-Triangularize once with inequalities and they seemed to slow it down so much I didn't wait for a solution, on a problem smaller than yours.

@MaPal93 In general, Maple doesn't handle assumptions particularly well, so I usually leave them out, or just add "assuming" clauses when I really need them. For solve itself, you need "useassumptions" as a solve option for the assuming clause to take effect - see the gamma file I uploaded for an example. For the polynomial systems if you look at RegularChains:-Triangularize that PolynomialSysyem is using, you can see how to more specifically add inequalities, but it might be faster to process them after as needed.

@MaPal93 Because a set cannot have duplicates, if you do a substitution that reduces nops from 3 to 2, then two of these equations became equal. (To see which are equal you could do the substitution on the individual equations and then compare them, by seeing if simplification of their differences was zero.)

I analysed the quadratic solution of the gamma equation, and found there are no all-positive solutions.

Edit: forgot the file:

220523_SOLUTION_reducedform_calibrations2.mw

@dharr the positive solution had lambda__1 = lambda__2 (some others did too). If there is some reason why this will be always true, then the system of 3 equations reduces to two, and can be solved extremely quickly with solve.

EqN2.mw

@MaPal93 I agree with your analysis, that the 5th solution has a unique positive solution, and that the first solution has an inifinite number of real solutions but none are positive. 

For the discriminant solving for when it is zero for a quadratic finds one real root, but positive discriminant gives other real roots. (In the general case it is zero if and only if a polynomial has roots of multiplicity greater than 1.)

Edit: note that solution 1 is just the common factor.

@davimon L^2 + L - (element in G) = 0 is still a polynomial with coefficients in G, so can be solved the same way.

Download GF24solve2.mw

@dharr 

Download GF24solve.mw

@davimon Sorry, I was working with your code outside the procedure. Outside the error message is that lambda is not a number in the field. Inside I don't really understand the error message. Usually I would use "uses packagename;"  in a procedure, but the GF package doesn't lend itself to this usage easily.

(See above edit, but I'm not sure yet how it applies to your case). But for sure solve is not going to work for your problem. [Edit -see below for how to this]

@davimon G has only numbers and does the field operations on them, it doesn't have variables. Solving lambda^2-(number in G) = 0 means solving a polynomial whose coefficient field is G. 

Edit: Factor and other commands can handle Galois Fields (but this isn't done through the G:-GF() construct.) From the Factor help page: "The call Factor(a, K) mod p and computes the factorization of the polynomial a over the Galois field (finite field) defined by the algebraic extension K over the integers mod p a prime integer. The extension K is specified by a RootOf and must be irreducible over the integers mod p."

p.s., the usage is

use G in 
c:=b^2
end use;

@MaPal93 

Let me be clear about my goal. I am looking for a (i) closed-form, (ii) real, and (iii) positive solution. Even just one solution satisfying all three constraints would be fantastic.

Now, it seems that sol[1] (2 equations in the lambdas) and sol[5] (3 equations in the lambdas) do not look weird and are the most promising for further analysis.

sol[2..5] all involve solutions of high degree polynomials, which in general do not have solutions in terms of radicals. So when you put your parameters in, you will not likely get symbolic solutions. So sol[1], which is surprisingly just a quadratic is your best chance. 

However, for the case with all parameters =1, an exhaustive analysis proves that there are no positive solutions for any of the 5 cases:
triade.mw

My questions:

1) How did you find the unique and real but negative lambdas from sol[5]?

In each case, you solve the third equation for lambda__1, substitute it into the second equation and solve it for lambda_2 and then substitute these solutions into the first equation and solve for lambda__3 (this is the backsubstitution on the triangularized system). See details in the file above. 

2) You showed me one example in which also sol[1] produces a real but negative solution. In the same comment you mentioned "but real and positive is harder". Why is it the case? If mathematically real and positive solutions are not ruled out a priori, is there a more systematic way to find them? (even just a loop which tries many substitutions and stops as soon as we find one would do, I guess?)

I did it in detail on the attached. With parameters it will be harder, but some things can help in this sort of analysis, Descartes rule of signs etc.

3) It's not clear to me the benefit of dividing out the common factor in the system with all parameters normalized to 1 (since we already found solutions and, factorized or not, these solutions would be the same). The benefit would be clear if such removable common factor existed in the uncalibrated equations EqN := ((numer@evala@:-Norm@numer)'tilde'@eval)(MyEqs). How to verify this possibility?

I agree. I assume factorizing before solving speeds up the solution process. Factorizing may not work with parameters but you can try the factor() command to see.

4) Surely a naive question, but would a solution found as in 2) solve also the uncalibrated system? That is, would simplify(eval(EqN,sol)) still give me 0 for EqN := ((numer@evala@:-Norm@numer)'tilde'@eval)(MyEqs) instead of EqN := ((numer@evala@:-Norm@numer)'tilde'@eval)(MyEqs,P='tilde'1)?

Seems very unlikely.

5) Related to 4). My end goal is to study how lambda_1, lambda_2, and lambda_3 vary with my parameters. Is it legit to pick a real and positive lambda_2 and lambda_3 (found as in 2)) and plug them back into uncalibrated EqN[1] (and solve for lambda_1), then pick lambda_1 and lambda_3 and plug them back into uncalibrated EqN[2] (and solve for lambda_2), finally, pick lambda_1 and lambda_2 and plug them back into uncalibrated EqN[3] (and solve for lambda_3)? See my original script 160523_stylized_triade.mw: the first equation is for lambda_1, second for lambda_2, and third for lambda_3.

I'm not sure I understand this, but in general the solution for all parameters=1 is not going to generalize to the case with parameters.

Really thank you for helping me out with this.

You're welcome. There are some interesting symmetries that might help you, For example EqN[2] and EqN[3] are the same if lambda__1 and lambda__2 are exchanged (and this was true for sol[1].) Is this expected?