Why does solve() after sym('x', 'positive') return a negative number?

I'm trying to solve a system of equations. I have defined 4 variables to be positive, but solve() returns a negative value. As far as I understand sym('x', 'positive') command, the result of solve should be an empty set if the only solution is negative... Am I missing something? Here is my code:

syms l_1 s_1 s_2 l_2 clear;

l_1 = sym('l_1', 'positive');

l_2 = sym('l_2', 'positive');

s_1 = sym('s_1', 'positive');

s_2 = sym('s_2', 'positive');

ft_1 =  10;
ft_2 = 10;
g_1 = 10;
beta = 0.9;
lambda = 0.6;
gama_0 = 0;
gama_1 = 0.5;
gama_2 = 0.2;
gama_3 = 0.2;
gama_4 = 0.2;
a = 1;
alfa_0 = 1;
alfa_1 = 0.5;
alfa_2 = 0.5;
L = 8;
delta_1 = 2;
delta_2 = 2;

V1 = log(ft_1 + g_1 + l_1) - delta_1*(l_1)^2 - delta_2*(s_1)^2+ beta*(log(ft_2 + lambda*(gama_0 + gama_1*s_1 + gama_2*l_1 + gama_3*a) + l_2)  - delta_1*(l_2)^2 - delta_2*(s_2)^2) + beta^2*(log(L*alfa_0*((0.5*(gama_0 + gama_1*s_1 + gama_2*l_1 + gama_3*a) + 0.5*(gama_0 + gama_1*s_2 + gama_2*(l_1 + l_2) + gama_3*a + gama_4*(gama_0 + gama_1*s_1 + gama_2*l_1 + gama_3*a)))^alfa_1)*(l_1 + l_2)^alfa_2)- delta_1*(L)^2);

pVps_1 = diff(V1,s_1)

pVpl_1 = diff(V1,l_1)

pVps_2 = diff(V1,s_2)

pVpl_2 = diff(V1,l_2)

[xl1, xl2, xs1, xs2] = solve(pVpl_1, pVpl_2,pVps_1, pVps_2, l_1, l_2, s_1, s_2) 

Results that I get are:

xl1 = -0.15635731765851100895734570889025

xl2 = -0.22318136524599875462767864215809 

xs1 = 0.13647789503117958882560294058743 

xs2 = 0.12006738235828457981821794454359

Answer 1

[Sean de Wolski]

I recommend against using ls and 1s (ELs and ONEs) as your variable names since it's very hard to read.

The solver is returning negative answers because it cannot find an analytical solution. Once it cannot find an analytical solution it uses a numeric solver that does not obey the constraints.

You can explicitly call the numeric solver using:

evalin(symengine, 'numeric::fsolve([pVpl_1, pVpl_2,pVps_1, pVps_2],[l_1=0..infinity,l_2=0..infinity,s_1=0..infinity,s_2=0..infinity])')

But this is failing for me.

[Tomaz]

But I'm solving for l_1 etc. which I defined to be positive... If I for example define x = sym('x', 'positive') and try to solve:

solve(x+4, x)

I get: ans =

[ empty sym ]

Thus, I would expect something similar above...

[Sean de Wolski]

see update