Optimisation in swi prolog

Say I want to find argmax(x,y,z) -1/2(20x^2+32xy +16y^2)+2x+2y.

subject to: x>=0, y>=0,z>=0 and -x-y+z =0.

I know the partial derivatives being set to 0 is :

-20x-16y+2=0 and -16x-16y+2 =0

so we could have x= 0 and y =1/8 and z=1/8.

How would I do this in Swi-prolog? I see that there is library simplex for linear solving, but this is a quadratic problem but the partial derivatives are not. (I am a bit confused!)

This is what I have:

:- use_module(library(simplex)).

my_constraints(S):-
 gen_state(S0),
 constraint([-20*x, -16*y] = 0, S0, S1),
 constraint([-16*x,-16*y] = 0, S1,S2),
 constraint([x] >= 0,S2,S3),
 constraint([y] >= 0,S3,S4),
 constraint([z] >= 0,S4,S5),
 constraint([-x-y+z] = 0,S5,S).

?- my_constraints(S), variable_value(S,x,Val1),variable_value(S,y,Val2).
false.

Solution

There are several issues here. First, just to get this out of the way: library(simplex) can only handle linear constraints. So yes, it cannot—at least not directly—be used to solve your actual problem.

But library(simplex) is often useful regardless, and so I want to quickly point out the following:

variable_value/3 only works on the solved tableau. This means that you must have invoked maximize/3 first.

For example:
```
?- my_constraints(S), maximize([x,y], S, Max), variable_value(Max, x, X).
S = ...,
Max = ...,
X = 0.
```
Note that you must change the final goal of my_constraint/1 to constraint([-1*x, -1*y,z] = 0, S5, S) to conform to the syntax required by this library.

That being said, let us now get to the core of the issue: There are well-known ways to iteratively solve quadratic optimization problems, using a series of linear programs and reasoning about gradients to get closer to a solution. Thus, library(simplex) can indirectly still be used to solve your problem.

In particular, check out the method of steepest ascent available from miscellaneous programs. It includes a small symbolic derivative calculator written in Prolog. Yes, it's "symbolic" ;-)

Plugging in your task, I get:

?- maximize(- 0.5*(20*x(1)^2 + 32*x(1)*x(2) + 16*x(2)^2) + 2*x(1) + 2*x(2),
   [[-1,0,0],
    [0,-1,0],
    [0,0,-1],
    [-1,-1,1],
    [1,1,-1]],
   [0,0,0,0,0],
   [0,0,0], Max).
Max = [4.298588509886033e-17, 0.125, 0.12500000000000006] ;
false.

Which is, up to the unbearable nastiness of floating point arithmetic, something that I hope you can work with.