How to format output from a GLPK Solver

Question

Hi I have a dataframe of 28 student names as indicies, and 12 book titles as columns. I am trying to run a pyomo solver to get the students into groups of 3 or fewer, and hope to give them their highest rated book choice.

My ask is how do I get the output to look like the initial dataframe rather than

('Student 1', 'Title1') :     0 :   1.0 :     1 : False : False : Binary
('Student 1', "Title2") :     0 :   0.0 :     1 : False : False : Binary
('Student 1', 'Title3') :     0 :   0.0 :     1 : False : False : Binary
('Student 2', 'Title1') :     0 :   0.0 :     1 : False : False : Binary
('Student 2', 'Title2') :     0 :   1.0 :     1 : False : False : Binary
('Student 2', 'Title3') :     0 :   0.0 :     1 : False : False : Binary
('Student 3', 'Title1') :     0 :   1.0 :     1 : False : False : Binary

etc.

Desired Output:

Name       Title1   Title2   Title3
Student 1     1       0        0
Student 2     0       1        0
Student 3     1       0        0 

This is what I have so far

model = ConcreteModel()

#Variables
Rank = rankings.keys()
Name = rankings.index

model.x = Var(Name,Rank,domain=Binary)

#ObjFct

model.rank= Objective(expr=sum([model.x[name,rank]*rankings.loc[name,rank] for name in Name for rank in Rank]), sense=minimize)

#Constraints
model.cons=ConstraintList()

#OneBook
for name in Name:
     model.cons.add(sum([model.x[name,rank] for rank in Rank]) == 1)
    
#GroupSize
for rank in Rank:
    model.cons.add(sum([model.x[name,rank] for name in Name]) <= 3)

#Every student in a group
for name in Name:
    model.cons.add(sum([model.x[name,rank] for rank in Rank]) == 1)
    
SolverFactory('glpk').solve(model)

model.display()