Here is one possible solution:
import numpy as np
# The physical constants set to unity in the Planck unit system:
# c, G, hbar, 1/(4.pi.eps0), kB
constants = np.array([299792458, 6.67384e-11, 1.054571726e-34,
8.9875517873681764e9, 1.3806488e-23])
dimensions = ('length', 'mass', 'time', 'charge', 'temperature')
si_units = ('m', 'kg', 's', 'C', 'K')
M = np.array([[1, 0, -1, 0, 0],
[3, -1, -2, 0, 0],
[2, 1, -1, 0, 0],
[3, 1, -2, -2, 0],
[2, 1, -2, 0, -1]])
# Invert the matrix to find the Planck quantities in SI units
Minv = np.linalg.inv(M)
planck_set = np.prod(constants**Minv, axis=1)
for dimension, planck_quant, si_unit in zip(dimensions, planck_set, si_units):
print('Planck {} = {} {}'.format(dimension, planck_quant, si_unit))
Which gives:
Planck length = 1.6161992562119142e-35 m
Planck mass = 2.1765092531302005e-08 kg
Planck time = 5.391060425582535e-44 s
Planck charge = 1.8755459567942785e-18 C
Planck temperature = 1.4168331315099011e+32 K
(Note that not all of these digits are significant).