Stuart,
I think you can do this as follows setting solveVariable to p3[1] --> the variable you want to solve for.
eqns = {p2[1] == p1[1] + p1[2], p2[2] == p1[3] + p1[4];
p2[3] == p1[5] + p1[6], p2[4] == p1[7] + p1[8],
p3[1] == p2[1] + p2[2], p3[2] == p2[3] + p2[4],
p4[1] == p3[1] + p3[2]}
solveVariable = p3[1]
The code below grabs all of the variables in the equations and sets them equal to vars. (this will only run in Version 10 or 11, otherwise do a DeleteCases to remove the solveVariable)
vars = DeleteDuplicates[Cases[eqns, _[_], Infinity]] /.
solveVariable -> Nothing
{p2[1], p1[1], p1[2], p2[3], p1[5], p1[6], p2[4], p1[7], p1[8], p2[2],
p3[2], p4[1]}
Now create all possible sets of variables (excluding the solveVariable) and then add the solveVariable to every set. Only create sets up to the number of equations - 1 because you can't have more variables than you have equations.
solveVarsSets =
Map[Prepend[solveVariable], Subsets[vars, Length[eqns] - 1]]
Run the Reduce to get your equations and delete the empty solutions and all duplicates. Only keep the equations in the form of the solveVariable == ....
answers =
DeleteDuplicates[
DeleteCases[
Map[Cases[Reduce[eqns, #], Equal[solveVariable, _]] & ,
solveVarsSets], {}]]
to get all possible equations that you want:
{{p3[1] == -p3[2] + p4[1]}, {p3[1] == -p2[3] - p2[4] + p4[1]}, {p3[
1] == p2[1] + p2[2]}, {p3[1] ==
p1[1] + p1[2] + p2[2]}, {p3[1] == -p1[5] - p1[6] - p2[4] +
p4[1]}, {p3[1] == -p1[7] - p1[8] - p2[3] + p4[1]}, {p3[
1] == -p1[5] - p1[6] - p1[7] - p1[8] + p4[1]}}
I hope this helps
Regards
Neil