Finding the maximum of a multivariate polynomial of degree one I need to find the global maximum of the function
\begin{align}
f\left(x\right) & = p_1  \max\left(\sum a_{1i} x_{1i}, \sum b_{1i} x_{1i}\right) - \sum c_{1i} x_{1i} \\
&+\ldots \\
&+ p_n  \max\left(\sum a_{ni} x_{ni}, \sum b_{ni} x_{ni}\right) - \sum c_{ni} x_{ni}
\end{align}
where all the coefficients $a_{ji}, b_{ji}$ and $c_{ji}$ are nonnegative and each variable $x_{ji}$ lives in $[0,1]$.
The number of terms $n$ is between 500 and 100,000.
I have some flexibility in choosing this $n$.
The index $i$ runs over 30 terms, so each summation in $f$ contains 30 terms.
What are my options to solve this problem?
I think $f$ is a convex function, so the maximization is an integer programming problem. But, at the same time, it has a linear structure, so I wonder if it is better to view as a linear programming problem.  
I'm considering a branch-and-bound algorithm, but I don't see what will be good lower an upper bounding functions. 
A gradient descent may work, but I need to be sure that it founds the global maximum.
Edit
The variable $x_{ji}$ are not all independents.
Some variables are repeated. For example, 
\begin{align}
&x_{11}=x_{21}=x_{31}=x_{41} \\
&x_{51}=x_{61}=x_{71}=x_{81} \\
&\ldots
\end{align}
and
\begin{align}
&x_{21}=x_{22} \\
&x_{23}=x_{24} \\
&\ldots
\end{align}
 A: You are trying to maximize a convex function over a convex set (in your case a slice of a high-dimensional cube. This has been studied a lot, and you can see this stackexchange discussion for references (though there are new ones almost daily, since this is of crucial importance in compressed sensing, distance matrix reconstruction, etc).
A: This is not an answer, and my expertise in tackling these problems
is very small.  My commentary above seemed to help, however;
Nicolas has a refined version of the problem in any event.  I have
another opinion to render which may inspire someone else to
post a solution.
The version of the problem looks like optimization of a
real valued function on a (piece of a) high dimensional
manifold, with two features that might be used in
developing a quick method for finding a global optimum:
the first is that function is composition of an order preserving
operator (sum from 1 to n) of other functions of "low dimension", and
the second is that these other functions are relatively simple in terms of
multivariable optimization; were it not for the relations between the labeled
variables, you could tell the computer to solve all the low dimensional problems,
tot up the results on an adding machine, and call it a day.
However (at this writing), the relations between the slices (my
temporary name for the projections of the manifold into 
compact subsets of various lower dimensional spaces) also
seem to be order preserving, and suggest promise of finding
quick optima.  The suggestion below is to solve the problem on
a slice, then on a related but initially disjoint slice, then laminate
the two slices together and adjust the result.
From the example posted in the question, I take the subterm of
f that involves variables x_1i and work on that domain of those
variables and call that slice S1; I do something similar with 1
replaced by 2 and get slice S2.  Lets solve on the two slices
and call the optimal inputs xb1 and xb2, assuming no
dependencies between the variables for S1 and the variables
for S2.
Time to laminate.  To simplify matters let's take a simple relation
of x_11=x_21.  (If one had x_11=-x_21, that would not be an order preserving
relation in my view; one might use a different coordinate system or some
other artifice to "get everything headed in the same direction", but I will shy
away from such difficulties for now.) If the first coordinates of xb1 and xb2
agree, then an appropriate stitching of the solutions should yield a
solution xb12 to the sum of the two terms on S12, 
a lamination of slices S1 and S2 (verification needed).  Otherwise, one
needs to combine xb1 and xb2 in some way to produce xb12.
I propose for each slice Sj to compute the differential cost d_jc
between the solutions already computed and a new optimum
under the additional constraint x_j1=c.  I won't embarrass myself by
working out the algebra and getting it wrong: instead I will
(hopefully not) embarrass myself by asserting d_jc is affine in c, or
linear in abs(c - xbj_j1).  Now we have a one dimensional linear
minimization which answer c is applied to the two restricted problems
on the slices, producing new values for xb1 and xb2 which agree on
the variables x_j1, which can be stitched together to form xb12.
Now the programme is to slice the (domain of the) manifold into
Si, solve the simple problems on these slices, and then laminate them back together,
hopefully in a computationally cheap way.  This is likely highly unoriginal as 
idea.  I hope it proves to be more motivational than unoriginal.
