Python3, 276 bytes

S=lambda b:sum(sum(i)==len(i)for i in b)
def f(b):
 q,d,U=[([],b)],{S(b):[[]]},[b]
 for m,b in q:
  for i in range(len(b[0])):
   B=[j[:i]+[int(not j[i])]+j[i+1:]for j in b];d[s]=d.get(s:=S(B),[])+[m+[i+1]]
   if B not in U:q+=[(m+[i],B)];U+=[B]
 return min(d[max(d)],key=len)

Try it online!

Pyth, 28

f@eo+/QN/Qm!dN_osZ^U2lhQTUhQ

Takes input in the form of a nested list, e.g.

[[0,0,1,0,0],[0,0,1,0,0],[0,0,1,1,0],[0,0,0,1,1],[0,0,0,1,1]]

Gives output 0-indexed, e.g.

[2]

^U2lhQ: Generates all possible lists of 0s and 1s of the right length.

_osZ: Orders these lists from most 1s to least.

+/QN/Qm!dN: Counts how many times each list (N) and its inverse, 0s and 1s swapped (m!dN) occur in the input. The former corresponds to a series of flips leaving all zeros, the latter to leaving all ones.

eo: Orders the list by the above key, and takes its last element, which will be the result with the most matching columns, and among them the one with the least ones.

f@ ... TUhQ: Converts this list of 1s and 0s to a list of indices to be flipped.

For 1-indexing, change the d to a k, then put mhd at the beginning.

Python 2, 234

from itertools import *
A=input()
n=len(A[0])
R=range(n)
S=(0,)
for p in[q for i in R[:-1] for q in combinations(R,i)]:
    s=[sum([(l[q]+(q in p))%2 for q in R])for l in A]
    m=s.count(n)+s.count(0)
    if m>S[0]:S=(m,p)
print S[1]

Input is a list of lists:

[[0,0,1,0,0],[0,0,1,0,0],[0,0,1,1,0],[0,0,0,1,1],[0,0,0,1,1]]

Output is a tuple of flips using python indexing from 0:

(2,)

If the output must be indexed from 1, then the code is 272 characters with 0 denoting no flips:

print 0 if len(S[1])==0 else [p+1 for p in S[1]]

Mathematica - 76

{
{P, P, T, P, P},
{P, P, T, P, P},
{P, P, T, T, P},
{P, P, P, T, T},
{P, P, P, T, T}
}/.{P -> 1, T -> 0};
First@MaximalBy[
  Subsets@Range@Length[%],
  MapAt[1-#&,%,{All,#}]~Count~{1..}&
]

CJam, 37

q~_{:!}%+:T{T1$a-,\0-+}$0={U):U*}%0-`

Input format:

[[0 0 1 0 0] [0 0 1 0 0] [0 0 1 1 0] [0 0 0 1 1] [0 0 0 1 1]]

APL, 37

{x/1+⍳⍴x←y[↑⍋(+/∘.≠⍨2⊥⍉y),¨+/y←⍵⍪~⍵]}

Example:

{x/1+⍳⍴x←y[↑⍋(+/∘.≠⍨2⊥⍉y),¨+/y←⍵⍪~⍵]} 5 5⍴0 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 1 1

Tested here.

Haskell, 98

g s=snd$maximum[((sum[1|o<-s,o==r||o==map(1-)r],-sum r),[i|(i,1)<-zip[1..]r])|r<-s++map(map(1-))s]

input format is a list of list of ints. you can use a string version for testing:

gg = g . map (map (\c -> case c of 'T' -> 0 ; _ -> 1) ) . lines

_{oh no the spaces! it hurts!}

this works by iterating over the rows and calculating the number of points we will get if we flipped the columns in a way that this row gains us a point.

first thing to notice is that flipping the row to either all True or all False doesn't matter, because the grids will be exactly the inverse of each other and thus will have the exact same score.

the way we calculate the count when a given row gains a point is so: we iterate over the rows again, and sum the points each row gives us, using the fact that they do exactly when the rows are either identical or the exact inverse.

for example, if the row we are flipping is TPPTP and the current row we are iterating over is PTTPT or TPPTP then the row gains us a point, but when it's any other row, it doesn't gain us any points.

CJam, 53 51 bytes

l~z:X,,La\{1$f++}/{,}${X\{_X=:!t}/z{_&,(},,}$0=:)S*

This reads a two-dimensional array of 0s and 1s from STDIN. E.g. the example in the question would be

[[0 0 1 0 0] [0 0 1 0 0] [0 0 1 1 0] [0 0 0 1 1] [0 0 0 1 1]]

Test it here.

This first gets all possible subsets of columns, in order of increasing length, then performs the flips for each subsets and sorts them by how many rows still have both 0s and 1s in them. Finally, we just return the first such subset. This relies on the sorting being stable, such that the initial order of increasing length takes care of the tie-breaker.