% META example
% Add2:  FF + Gain + Integrator --> Saturation17 [20K,75K]
% Add1:  FLOW_SP - FLOW_FB --> Gain2 of Gain=10
% RST,N_init,Gain2 -> Integrator --> 
% N_init = 0
% Add1 ->  Gain2 -> Gain of gain=50
% SPEED_SET = Saturation of FF + 50 * 10 * (SP-FB) + inte;
% Final Model is:
% d/dt inte = e;
% d/dt e = - alpha * ( 500*e + 10*inte )
% Therefore, Kp = 500, Ki = 10, 
% Template: c * e^2  + d * inte^2
% Template for assumption:  d/dt fb = alpha * control,
% where alpha varies within some interval [aa,bb];

% aa = bb > 0 is the only solution; i.e. we don't have
% any Lyapunov function that works for more than ONE dynamics!!
% Even if we generalize the form of the Lyapunov function a bit, 
% we still get aa = bb as the only solution.
% Answer obtained using "slfq -N 100 pid-ag.slfq > pid-ag.slfq.out"
% We also tried adding precondition of V being bounded by 1.
% Even replacing > by < in the precondition of the implication
% has the same answer, though the answer is difficult to obtain:
% using slfq on first 15 facts, we can deduce d = 10 c aa = 10 c bb
% Changing  (c*e^2 + d*inte^2 < a) by a := 1

% Changing problem: template is (c*e + d*inte)^2.  
% There was a mistake in the old version, hence results were wrong.
% Under the assumption that (c*e + d*inte)^2 = 1 
% we get aa > 0 and bb > aa.;  without the assumption, it was FALSE
% But value of c and d found was c = d = 0 !
% For a particular alpha, there is one c for every d, hence no 
% hope for V that works for more than one alpha.

% If I use template: (-50*d*e + d*intwe)^2 + b*intwe^2
% and use forall alpha > 0, (and no aa,bb), then
% then virtual+qe+slfq yields condition on b,c,d as:
% [ c = 0 /\ 2500 b + c^2 - 100 d c + 2500 d^2 <= 0 ]
% which simplifies to c=0 and 2500 b + 2500 d^2 <= 0 ...b + d^2 <= 0
% which does not give a good solution.
% Changing to b*e^2 also gives no nontrivial solution.
% Restricting alpha to 1/1000 to 1/100 also does not help.

% Finally, the interesting case is when I assume
% that the integrator output is saturated to +-1.
% SO, now the system becomes hybrid (switched) depending
% on whether inte has reached its saturation value.
% We find a COMMON Lyapunov function for each mode
% of the form c*e^2 + d*inte^2.
% We get the condition d <= 510 * aa* c 
% So, there is a V that works for all alpha !!

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
rlqepcadl(200*10^3);  % The size of the prime list used by QEPCAD. The
on rlverbose;
on time;
%off rlqefb;		% No fallback QE

% phi1 := ( c*e^2 + d*inte^2 + b*e*inte >= 0 );
% (c*e+d*inte)^2 -->

phi2 := 
  ( e^2 >= 1 and (inte^2 < 1 or (inte^2 = 1 and e*inte < 0)) and 
	aa <= alpha and alpha <= bb )  impl
  ( ( -2*e*c*alpha*(500*e+10*inte) + 2*d*inte*e < 0 ) );

phi3 := 
  ( ((e >= 1 and inte = 1) or (e <= -1 and inte = -1)) and 
	aa <= alpha and alpha <= bb )  impl
  ( ( -2*e*c*alpha*(500*e+10*inte) + 2*d*inte*0 < 0 ) );


% phi2 := 
  % ( ( e > 0 and inte > 0 and aa <= alpha and alpha <= bb) ) impl 	% This works...
  % ( ( -2*e*c*alpha*(500*e+10*inte) + 2*d*inte*e < 0 ) );
  
% phi3 := 
  % ( ( e < 0 and inte > 0 and aa <= alpha and alpha <= bb) ) impl 	% This works...
  % ( ( -2*e*c*alpha*(500*e+10*inte) + 2*d*inte*e < 0 ) );
% 
% phi4 := 
  % ( ( e < 0 and inte < 0 and aa <= alpha and alpha <= bb) ) impl 	% This works...
  % ( ( -2*e*c*alpha*(500*e+10*inte) + 2*d*inte*e < 0 ) );
% 
% phi5 := 
  % ( ( e > 0 and inte < 0 and aa <= alpha and alpha <= bb) ) impl 	% This works...
  % ( ( -2*e*c*alpha*(500*e+10*inte) + 2*d*inte*e < 0 ) );
% 
phi := all({alpha,e,inte}, c > 0 and d > 0 and aa > 0 and bb >= aa and phi2 and phi3 );
% phi := all({alpha,e,inte}, aa > 0 and aa < bb and phi2 ) ;

% psi := rlqepcad( phi, "pid-ag2.qepcad" )$
% psi := rlqepcad( phi );

psi := rlqe( phi )$
psi2 := rlqepcad( rlex(psi), "pid-ag2.qepcad" )$

off nat;	% Output same syntax as input (for reduce)
out "pid-ag2.out";
psi;
psi2;
shut "pid-ag2.out";
 
end;