In practice (at least in the rates world), $\beta$ is preset and $\alpha$ is solved for to calibrate to the atm vols $\sigma_{ATM}$ (which are the most liquid and reliable of the market data available). For instance, in the case for normal vols and assuming a normal distribution of the fwds, $\beta=0$ then $$σ_{N,ATM}=α\left(1+\frac{2−3ρ^2}{24}ν^2T\right).$$ The $\rho,\nu$ parameters are then obtained via the sort of optimization routines you describe in order to incorporate the skew.

Edit: I have added a simple SABR calibration routine (employing RSS) I use to illustrate what I mean below.

    def sabr_calibration(swvolcube,a,b,calculation_date,alphas,tolerance):
    print('for '+str(a.normalized())+str(b.normalized())+' optimizing alpha to fit sabr to atm vols ...')
    T=SABR(swvolcube,a,b,calculation_date,'T')
    f=SABR(swvolcube,a,b,calculation_date,'f')
    atm_sabr_vol=SABR(swvolcube,a,b,calculation_date,'atm_sabr_vol')
    atm_vol=SABR(swvolcube,a,b,calculation_date,'atm_vol')
    atm_smile=SABR(swvolcube,a,b,calculation_date,'atm_smile')
    strike_spreads=SABR(swvolcube,a,b,calculation_date,'strike_spreads')
    atm_strike_set=[f+i for i in strike_spreads]
    beta=SABR(swvolcube,a,b,calculation_date,'beta')
    nu=SABR(swvolcube,a,b,calculation_date,'nu')
    rho=SABR(swvolcube,a,b,calculation_date,'rho')
    cubic0=-atm_vol*(f**(-beta))
    cubic1=(1+((2-3*(rho**2))/24)*(nu**2)*T)
    cubic2=(rho*beta*nu*T)/(4*(f**(1-beta)))
    cubic3=(beta*(beta-2)*T)/(24*(f**(2-2*beta)))
    coeff=[cubic3,cubic2,cubic1,cubic0]
    roots=[np.roots(coeff)[i] for i in range(0,len(np.roots(coeff)))]
    positive_roots=[i for i in roots if i>0]
    positive_roots.sort()
    root=positive_roots[0]
    dummy_alpha = alphas.loc[str(a.normalized()).lower(), str(b.normalized()).lower()]
    dummy_nu = nus.loc[str(a.normalized()).lower(), str(b.normalized()).lower()]
    dummy_rho = rhos.loc[str(a.normalized()).lower(), str(b.normalized()).lower()]
    if abs(atm_vol-atm_sabr_vol)>tolerance:
        dummy_alpha.setValue(root)
        sabr_calibration(swvolcube, a, b, calculation_date,alphas,tolerance)
    else:
        print('atm vol = ' + str(atm_vol))
        print('atm sabr vol = '+str(atm_sabr_vol))
        print('optimized alpha = '+str(dummy_alpha.value()))
        print('error is = ' + str(abs(atm_vol - atm_sabr_vol)))
        print('calibrating nu and rho ...')
        params = np.array([nu, rho])
        def calib(params):
            vols = np.array([
                ql.sabrVolatility(strike, f, T, dummy_alpha.value(), beta, *params, vol_type)
                for strike in atm_strike_set
            ])
            return ((vols - np.array(atm_smile)) ** 2).mean() ** .5
        cons = (
            {'type': 'ineq', 'fun': lambda x: 0.999 + x[1]},
            {'type': 'ineq', 'fun': lambda x: 0.999 - x[1]},
            {'type': 'ineq', 'fun': lambda x: x[0] - 1e-15},
        )
        result = minimize(calib, params,constraints=cons)
        new_params = result['x']
        nu = new_params[0]
        rho = new_params[1]
        dummy_nu.setValue(nu)
        dummy_rho.setValue(rho)
        atm_sabr_vol = SABR(swvolcube, a, b, calculation_date, 'atm_sabr_vol')
        atm_vol = SABR(swvolcube, a, b, calculation_date, 'atm_vol')
        if abs(atm_vol - atm_sabr_vol) > tolerance:
            print('re-calibrating ...')
            sabr_calibration(swvolcube, a, b, calculation_date, alphas, tolerance)

    return dummy_alpha.value(),nu,rho,beta