Intended behavior for dimension errors on `SymbolicDimensions`

I would like to hear the opinions of @odow @trulsf for JuMP and maybe @ChrisRackauckas or @YingboMa for ModelingToolkit.jl.

What would you prefer the behavior to be when you try to add two quantities that have incompatible `SymbolicDimensions`?

First, recall that SymbolicDimensions is a way to avoid automatically converting to SI base units:

For example:

```julia
julia> q = 4us"kPa * Constants.c"
4.0 kPa c

julia> sqrt(q)
2.0 kPa¹ᐟ² c¹ᐟ²
```

Only when you call `expand_units` would it convert from `SymbolicDimensions` to `Dimensions`:

```julia
julia> expand_units(q)
1.199169832e12 kg s⁻³
```

Right now when you try to add two incompatible quantities, the behavior is:

```julia
julia> 1us"km/s" + 1us"m/s"
ERROR: DimensionError: 1.0 s⁻¹ km and 1.0 m s⁻¹ have incompatible dimensions
```

whereas with regular `Dimensions`, you get:

```julia
julia> 1u"km/s" + 1u"m/s"
1001.0 m s⁻¹
```

**Questions for you:** is this behavior good as-is? Or should there be some way to automatically convert to SI units when this happens, but stay in `SymbolicDimensions` (to avoid type instability)?

---

Another option would be to perform all computations with `Dimensions`, and then display by converting back to `SymbolicDimensions`. For example:

```julia
julia> output_units = us"km/hr";

julia> a = 15us"km/hr"; b = 30us"m/s";

julia> raw_output = expand_units(a) + expand_units(b)
34.166666666666664 m s⁻¹

julia> as_units(q, units) = Quantity(ustrip(q / expand_units(units)) * ustrip(units), dimension(units));

julia> as_units(raw_output, output_units)
122.99999999999999 km hr⁻¹
```

Start with the SymbolicDimensions arithmetic and DimensionError behavior shown in the examples, then compare the regular Dimensions and expand_units paths. The issue names no files or tests; done means reaching agreement on incompatible additions and validating the chosen behavior, including the type-stability tradeoff.

Written by the indexing model from the issue text.