% Inverted Pendulum Example from EMSOFT 2010 paper. %% Dynamics: 3 MODES: +-0 % pos = vel % vel = +-2 - 15/100 omega omega theta + 5 theta - 15/2 theta theta theta -+ 2 theta theta + 15/100 omega omega theta theta theta % theta = omega % omega = 50 theta +- 20/3 - 1/2 omega omega theta - 100/3 theta theta theta -+ 10 theta theta + 3/4 omega omega theta theta theta +- 10/3 theta^4 %% Safe: 20*theta^2 <= 1 %% Template: [-theta^2 - b*omega^2 + c >= 0] %% before forall omega can be elimiated, we already see %% bounds on c in the formula; we pick c = 1/100 arbitrarily within %% the bounds and then eliminate omega to get condition on b as: %% 4801 b - 300 >= 0 on echo; load_package redlog; rlset r; load_package qepcad; on rlqefbqepcad; % use QEPCAD as the fallback QE rlqepcadn(100*10^6); % The number of cells used by QEPCAD. The % default is (only) 1 * 10^6. rlqepcadl(200*10^3); % The size of the prime list used by QEPCAD. The % default is (only) 2 * 10^3. on rlverbose; on time; off rlqefb; % No fallback QE phi1 := (( 20*theta >= -1 and 1 >= 20*theta and omega=0 ) impl (-1*theta^2 - omega^2*b + (1/100) >= 0)); phi2 := ( (-1*theta^2 - omega^2*b + (1/100) >= 0) impl ( 1 >= 20*theta^2 ) ); phi3 := ( (-1*theta^2 - omega^2*b + (1/100) = 0 ) impl ( -3/2*theta^3*omega^3*b + 200/3*theta^3*omega*b + theta*omega^3*b - 100*theta*omega*b - 2*theta*omega >= 0 OR 20/3*theta^4*omega*b - 3/2*theta^3*omega^3*b + 200/3*theta^3*omega*b - 20*theta^2*omega*b + theta*omega^3*b - 100*theta*omega*b - 2*theta*omega + 40/3*omega*b >= 0 OR -20/3*theta^4*omega*b - 3/2*theta^3*omega^3*b + 200/3*theta^3*omega*b + 20*theta^2*omega*b + theta*omega^3*b - 100*theta*omega*b - 2*theta*omega - 40/3*omega*b >= 0 ) ) ; phi := phi1 and phi2 and phi3; psi := rlgqe( all({pos, vel, theta, omega}, phi) ); psi := second psi; psi2 := rlqepcad( rlex(psi), "inverted-pend.qepcad" ); off nat; % Output same syntax as input (for reduce) out "inverted-pend.out"; psi; psi2; shut "inverted-pend.out"; end;