How to convert quadratic to linear program?

I have an optimization problem that has in the objective function 2 multiplied variables, making the model quadratic.

I am currently using zimpl, to parse the model, and glpk to solve it. As they don't support quadratic programming, I would need to convert this to an MILP.

. The first variable is real, in range [0, 1], the second one is real, from range 0 to inf. This one could without a problem be integer.

The critical part in the objective function looks like this:

max ... + var1 * var2 + ...

I had similar problems in the constraints, but they were easily solvable.

How could I solve this kind of a problem in the objective function?

Models in this form are actually called bilinear optimization problems. The typical approach to linearizing bilinear terms is through something called the McCormick envelope.

Consider variables x and y, where you want x*y in the objective of your maximization problem. If we assume x and y are bounded by xL <= x <= xU and yL <= y <= yU, then we can replace x*y with w, an upper bound for the quantity, with the following constraints (you can see the derivation here):

w <= xU*y + x*yL - xU*yL
w <= x*yU + xL*y - xL*yU
xL <= x <= xU
yL <= y <= yL

Note that these constraints are all linear in the decision variables. There are corresponding lower bounds in the McCormick envelope, but since you're maximizing they're unimportant in your case.

If you want a tighter bound on x*y, you can split up the interval on one of the variables (I'll use x here) into ranges [xL1, xU1], [xL2, xU2], ..., [xLn, xUn], introducing auxiliary continuous variables {x1, x2, ..., xn} and {w1, w2, ..., wn} as well as auxiliary binary variables {z1, z2, ..., zn}, which will indicate which range of x values was selected. The constraints above could be replaced by (I'll show the index 1 case, but you would need these for all n indices):

w1 <= xU1*y + x1*yL - xU1*yL*z1
w1 <= x1*yU + xL1*y - xL1*yU*z1
xL*z1 <= x1 <= xU*z1

Basically you will have x1=0 and w1 <= 0 whenever z1=0 (aka this part of the range is not selected), and you will have the normal McCormick envelope if z1=1 (aka this part of the range is selected).

The last step is to generate x and w out of range-specific versions of these variables. This can be done with:

x = x1 + x2 + ... + xn
w = w1 + w2 + ... + wn

The larger you make n, the more accurate of an estimation you will have for the bilinear term. However large values of n will affect the tractability of solving your model.

One final note -- you indicate that one of your variables in unbounded, but the McCormick envelope requires bounds on both variables. You should fix bounds, solve, and if your optimal value is at the boundary then you should re-solve with a different boundary.