Solution extended from answer @smichr where I rewrote the list comprehension with 'for' loops so I can immediately solve the imaginary part, and put in the result through subs to get the final_result and update the dict. Then I check if the dicts are unique before returning the result.

result = sp.solve(AllEquestions, AllVariables, dict=True)
want = set(AllVariables)
unk = []
for s in result: #loop through all solutions
    for key, value in s.items(): #loop through individual values
        if value.is_real is None and value.free_symbols&want: #check if value is img and is free_symbol e.g. value = 1.03077640640442*I*(10.0 - z4) + 16.25
            result_img_part = sp.solve(sp.im(value), want) #solve only the img part e.g. result_img_part = [{z4: I*im(z4) + 10.0, re(z4): 10.0000000000000}]
            for result_real_part in result_img_part[0].values(): #search for the real part in the solution
                if result_real_part.is_real:
                    final_result = value.subs(list(value.free_symbols)[0], result_real_part) #solve the function (value) to get the real value
                    s[key] = final_result # update the dictionary
    unk.append(s) 

#return unique solutions
res = []
for _dict in unk:
    if _dict not in res:
        res.append(_dict)

if result == []:
    raise Exception("No solution found!")
else:
    return res

Final thoughts. Code should be further improved, in the case there are multiple free_symbols to solve the imaginary part.