integralpwa.m

%% Look at function
ti = 0:0.001:8;
fi = min(min(2*ti,4),16-3*ti);
l = plot(t,fi);
grid on;hold on
l = plot(ti,0.001*cumsum(fi));

%% Manually derived
yalmip('clear')
sdpvar x
sdpvar f
region = binvar(3,1);

R1 = [0 <= x <= 2];
R2 = [2 <= x <= 4];
R3 = [4 <= x <= 10];
Model = [implies(region(1), [R1, f == x^2])
         implies(region(2), [R2, f == -4 + 4*x])
         implies(region(3), [R3, f == -28 + 16*x - 1.5*x^2])
         sum(region) == 1
         [0 <= x <= 10, -100 <= f <= 100]];

optimize(Model, -f)

%% Avoid quadratic equalities
sdpvar x1 x2 x3
Model = [implies(region(1), [R1, x1 == x, x2 == 0, x3 == 0])
         implies(region(2), [R2, x2 == x, x1 == 0, x3 == 0])
         implies(region(3), [R3, x3 == x, x1 == 0, x2 == 0])
         sum(region) == 1
         [0 <= [x x1 x2 x3] <= 10]];
f = x1^2 + (-4*region(2)+4*x2) + (-28+16*x3 - 1.5*x3^2);
optimize(Model, -f)

%% Lazy
sdpvar f
sdpvar z
f1 = 2*z; 
f2 = 4;   
f3 = 16-3*z;
q1 = int(f1,z,0,x);
q2 = int(f1,z,0,2) + int(f2,z,2,x);
q3 = int(f1,z,0,2) + int(f2,z,2,4) + int(f3,z,4,x);
Model = [implies(region(1), [R1, f == q1])
         implies(region(2), [R2, f == q2])
         implies(region(3), [R3, f == q3])
         sum(region) == 1
         [0 <= x <= 10, -100 <= f <= 100]];
optimize(Model, -f)

%% Trust but verify!
sdisplay(q1)
sdisplay(q2)
sdisplay(q3)


%% Even more lazy
f = int(f1,z,0,min(x,2)) + int(f2,z,2,min(x,4)) + int(f3,z,4,x);
optimize([0 <= x <= 10], -f)


%% Crazy (MILP approximation based on numerically computed integral)
f = interp1(ti,0.001*cumsum(min(min(2*ti,4),16-3*ti)),x,'sos2');
optimize([0 <= x <= 10],-f)

%% Just bad (nonlinear programming with on-the-fly computed integral)
sdpvar x
f = blackbox(@(x)(integral(@(z)(min(min(2*z,4),16-3*z)),0,x)),x);
optimize([0 <= x <= 10],-f);