Here's an attempt at automatic unit simplification. I'll give examples first, and implementation later.

First, define some target units:

SIbase = {"Seconds", "Meters", "Kilograms", "Amperes", "Kelvins", "Moles", 
          "Candelas", "Steradians", "Radians"};
SIderived = {"Newtons", "Pascals", "Joules", "Watts", "Coulombs", "Volts", 
             "Farads", "Ohms", "Webers", "Teslas", "Henries", "Lumens", "Lux"};
physicalconstants = {"SpeedOfLight", "PlanckConstant", "BoltzmannConstant", 
                     "ElementaryCharge", "ElectricConstant", "MagneticConstant", 
                     "GravitationalConstant", "JosephsonConstant"};

Express $\hbar^2/(m e)$ in terms of SI units:

Z = Quantity["ReducedPlanckConstant"]^2/(Quantity["ElectronMass"]*Quantity["ElectronCharge"]);
unitSimplify[Z, Join[SIbase, SIderived]]
(*    {7.61996423*10^-20 m^2 V}    *)

Express $h/k_B^2$ in SI units: who would have thought that the most concise way of expressing the unit is squared Kelvins per Watt?

unitSimplify[Quantity["PlanckConstant"]/Quantity["BoltzmannConstant"]^2,
             Join[SIbase, SIderived]]
(*    {6626070150000000000000000/1906191661201 K^2/W}    *)

Express the speed of light in funny units: there are four equivalently simple solutions,

unitSimplify[Quantity["SpeedOfLight"],
             {"Meters", "Seconds", "Hertz", "Wavenumbers"}]
(*    {299792458 m/s,
       299792458 m Hz,
       29979245800 per wavenumber per second, 
       29979245800 Hz/wavenumber}    *)

Express a meter in terms of physical constants:

unitSimplify[Quantity[1, "Meters"], physicalconstants]
(*    {2.4683*10^34 Sqrt[G] Sqrt[h]/c^(3/2)}    *)

implementation

First, a helper function: find the sparsest vector $\vec{x}$ of the underdetermined linear system of equations $\vec{x}\cdot A=\vec{b}$: (this is an NP-complete problem and I'm still looking for a more efficient implementation)

sparseSolve[A_, b_, n_] := 
  Quiet@Check[{#, LinearSolve[Transpose[A[[#]]], b]}, Nothing] & /@ 
    Subsets[Range[Length[A]], {n}]
sparseSolve[A_, b_] := Module[{n, solutions, shortsolutions},
  (* find the solutions with smallest number n of nonzero entries *)
  n = 1;
  While[(solutions = sparseSolve[A, b, n]) == {}, n++];
  (* of these, take the solutions that are smallest by 1-norm *)
  shortsolutions = MinimalBy[solutions, Norm[#[[2]], 1] &];
  (* convert to solution vectors *)
  SparseArray[Thread[Rule @@ #], Length[A]] & /@ shortsolutions]

The automated unit simplifier: simplify the units of Q by forming combinations of the units given in the list U. A list of equivalently simple solutions is returned:

unitSimplify[Q_Quantity, U_List] := 
  Module[{Uunitdimensions, unitdimensionslist, Uunitexponents, Qunitexponents,
          simplestunits},
    (* find the unit dimensions of each unit given in the list U *)
    Uunitdimensions = UnitDimensions[UnitConvert[#]] & /@ U;
    (* make a list of all the unit dimensions used here *)
    unitdimensionslist = Union @@ Uunitdimensions[[All, All, 1]];
    (* for each unit in U, make a list of exponents of the unit dimensions in *)
    (* the order of unitdimensionslist                                        *)
    Uunitexponents = Lookup[Rule @@@ # & /@ Uunitdimensions, unitdimensionslist, 0];
    (* for the desired unit Q, make a list of exponents of the unit dimensions *)
    (* in the order of unitdimensionslist                                      *)
    Qunitexponents = Lookup[Rule @@@ UnitDimensions[Q], unitdimensionslist, 0];
    (* find the simplest possible ways of expressing Qunitexponents as a linear *)
    (* combination of the rows of Uunitexponents                                *)
    simplestunits = Times @@ (U^#) & /@ sparseSolve[Uunitexponents, Qunitexponents];
    (* for each solution, convert Q to this unit *)
    UnitConvert[Q, #] & /@ simplestunits]