We present a method for fully automated generation of high quality robust proton treatment plans using hierarchical optimization. To fill the gap between the two common extreme robust optimization approaches, that is, stochastic and worst-case, a robust optimization approach based on the p-norm function is used whereby a single parameter, , can be used to control the level of robustness in an intuitive way. A fully automated approach to treatment planning using Expedited Constrained Hierarchical Optimization (ECHO) is implemented in our clinic for photon treatments. ECHO strictly enforces critical (inviolable) clinical criteria as hard constraints and improves the desirable clinical criteria sequentially, as much as is feasible. We extend our in-house developed ECHO codes for proton therapy and integrate it with a new approach for robust optimization. Multiple scenarios accounting for both setup and range uncertainties are included (13scenarios), and the maximum/mean/dose-volume constraints on organs-at-risk (OARs) and target are fulfilled in all scenarios. We combine the objective functions of the individual scenarios using the p-norm function. The p-norm with a parameter or result in the stochastic or the worst-case approach, respectively; an intermediate robustness level is obtained by employing -values in-between. While the worst-case approach only focuses on the worst-case scenario(s), the p-norm approach with a large p value ( ) resembles the worst-case approach without completely neglecting other scenarios. The proposed approach is evaluated on three head-and-neck (HN) patients and one water phantom with different parameters, . The results are compared against the stochastic approach (p-norm approach with ) and the worst-case approach, as well as the nonrobust approach (optimized solely on the nominal scenario). The proposed algorithm successfully generates automated robust proton plans on all cases. As opposed to the nonrobust plans, the robust plans have narrower dose volume histogram (DVH) bands across all 13 scenarios, and meet all hard constraints (i.e., maximum/mean/dose-volume constraints) on OARs and the target for all scenarios. The spread in the objective function values is largest for the stochastic approach ( ) and decreases with increasing toward the worst-case approach. Compared to the worst-case approach, the p-norm approach results in DVH bands for clinical target volume (CTV) which are closer to the prescription dose at a negligible cost in the DVH for the worst scenario, thereby improving the overall plan quality. On average, going from the worst-case approach to the p-norm approach with , the median objective function value across all the scenarios is improved by 15% while the objective function value for the worst scenario is only degraded by 3%. An automated treatment planning approach for proton therapy is developed, including robustness, dose-volume constraints, and the ability to control the robustness level using the p-norm parameter , to fit the priorities deemed most important.