Object Oriented Modeling with Objectica

Stephan Leibbrandt

({{, Welcome}})

({{Characteristics,  }, {Concepts of Ojectica,  }, {Heat Flow , Finite Differences}, {Traffic Dynamics, System of PDEs}, {Analog Computer Circuits , Coupling}, {Further Potentials,  }})

({{Characteristics, abstraction layer}, { , OO libraries}})

({{, ({{integrated}, {simple syntax}, {fast}, {code checks}})}})

({{Concepts of Ojectica, one class template -> many objects}, { , encapsulation of functionality}})

Standard Package

Get["Class`Class`"]

Class, Object and Access

Class[Person]:={
name=0,
Virtual.kind[]:="someone",
Person[name_]:=This.name=name,
info[]:=kind[]<>" called "<>name
}

per=New.Person["Miller"]

per.info[]

Inheritance and Polymorphism

Class[Employee,Person]:={
Private.salary=0,
Employee[name_,sal_]:=(Super[name];salary=sal),
Virtual.Override.kind[]:="worker with salary "<>ToString[salary]
};

emp=New.Employee["Smith",10000]

emp.info[]

Some Details

Clone

per2=CloneObject[per]

Save for External Use

Export[ToFileName[$TemporaryPrefix,"per"],Person,LicenseID->"L1234-5678"]

Export[ToFileName[$TemporaryPrefix,"per"],per]

Thread

{per,per2}.name

Variable Generation

Class[Var]:={
get[]:=$.x
};
Table[New.Var[],{5}].get[]

Straightforward Operator Usage

Class[Field]:={conts=0};
obj=New.Field[]
obj.conts++
obj.conts

Class[Field]:={conts={1,2,3}};
obj=New.Field[]
obj.conts[[2]]=27
obj.conts

Documentation

Class[Var]:={
generated=DateList[],
"Generate object symbol.".get[]:=$.x
};
Var

New.Var[]

({{Heat Flow, stepwise refinement -> iterative model development}, { , black box models -> nonexpert can build models}})

Idea

1D heat flow in layered temperature dependent medium

Graphics3D[{{Hue[1/3],Table[Cylinder[{{0,slice-1,0},{0,(1-1)*5+slice,0}},1],{slice,2}]},{Hue[2/3],Table[Cylinder[{{0,2+slice-1,0},{0,2+slice,0}},1],{slice,6}]},{Hue[1/3],Table[Cylinder[{{0,8+slice-1,0},{0,8+slice,0}},1],{slice,2}]}},BoxRatios->{1,10,1},AspectRatio->1/2,ViewPoint->{1.4,-2.9,1.},ViewVertical->{-1.,-0.4,9}]

({{Class, Object}, {({{Piece:, 0==(∂_x) (λ[θ[x]](∂_x) θ[x])==(∂_θ) λ[θ[x]]((((∂_x) θ[x]))^2)+λ[θ[x]](∂_(x,x)) θ[x]}, { , (∂_x) f[x+((Δx)/2)]≈((f[x+Δx]-f[x])/(Δx))}, { , (∂_(x,x)) f[x]≈(1/((Δx)^2))(f[x+Δx]+f[x-Δx]-2f[x])}, { , 0==(((∂_θ) λ[θ[x]](θ[x+Δx]-θ[x-Δx]))/(4((Δx)^2)))+λ[θ[x]](((θ[x+Δx]+θ[x-Δx]-2θ[x]))/((Δx)^2))}, { , 0==(∂_θ) λ[θ[x]](θ[x+Δx]-θ[x-Δx])+4λ[θ[x]](θ[x+Δx]+θ[x-Δx]-2θ[x])}}), create pieces}, {({{Materials:, ConstantConductivity}, { , LinearConductivity}, { , LDPE/HDPE}}), extract equations}, { , solve}})

Piece (Slices)

Class[Piece,Import:>{t,T0,T1}] := {
min=0,
max=0,
λ=0&,
Piece[min_Integer,max_Integer,λ_Function]:=(This.min=min+1;This.max=max;This.λ=λ),
Equations[]:=Table[0==D[λ[t],t]($.θ.(i+1)-$.θ.(i-1))+4λ[t]($.θ.(i+1)+$.θ.(i-1)-2$.θ.i)/.{t->$.θ.i,$.θ.0->T0,$.θ.11->T1},{i,min,max}]
};

Materials

λ == Constant

Class[ConstantConductivity,Piece] := {
ConstantConductivity[min_Integer,max_Integer,λ_Real]:=Super[min,max,λ&]
};

Class[MeanLDPE,ConstantConductivity] := {
MeanLDPE[min_Integer,max_Integer]:=Super[min,max,0.3125]
};

Class[MeanHDPE,ConstantConductivity] := {
MeanHDPE[min_Integer,max_Integer]:=Super[min,max,0.39375]
};

λ  ~ θ

Class[LinearConductivity,Piece] := {
LinearConductivity[min_Integer,max_Integer,basecond_Real,slopecond_Real]:=Super[min,max,basecond+slopecond(#-273.15)&]
};

Class[LDPE,LinearConductivity] := {
LDPE[min_Integer,max_Integer]:=Super[min,max,0.36875,-0.001125]
};

Class[HDPE,LinearConductivity] := {
HDPE[min_Integer,max_Integer]:=Super[min,max,0.44375,-0.001]
};

Layered Medium

rod=	  {New.HDPE[0,2],         New.LDPE[2,8],          New.HDPE[8,10]}
constrod={New.MeanHDPE[0,2],New.MeanLDPE[2,8],New.MeanHDPE[8,10]}

equations=Flatten[rod.Equations[]]/.{T0->293.15,T1->353.15}
constequations=Flatten[constrod.Equations[]]/.{T0->293.15,T1->353.15}

constraints=Thread[293.15≤Table[$.θ.i,{i,10}]≤353.15]

start=Thread[{Table[$.θ.i,{i,10}],Take[Table[i,{i,293.15,353.15,(353.15-293.15)/11}],{2,-2}]}]

Temperature Distribution

sol=FindMinimum[Prepend[constraints, #.#&[Last/@equations]],start]

constsol=FindMinimum[Prepend[constraints, #.#&[Last/@constequations]],start]

ListPlot[(Table[$.θ.i,{i,10}]/.Last[constsol])-(Table[$.θ.i,{i,10}]/.Last[sol])/.{T0->293.15,T1->353.15},Joined->True,PlotRange->All,Frame->True,Axes->False]

({{Traffic Dynamics, inheritance -> stepwise adding of functionality}, { , easy to change -> refinement driven by simulation results}})

Idea

distance of cars with finite reaction time

http://www.strassen.nrw.de/projekte/a2/rhynern.html

({{Class, Object}, {({{equations:, (Subscript[m, 0])*(Subscript[x, 0])''[t]==(Subscript[x, 0])[t], -> Car}, { , (Subscript[m, i])*(Subscript[x, i])''[t]==(Subscript[breakforce, i]),  }}), create cars}, {reaction time -> DelayedCar, ({{(Subscript[x, 0])[t]}, {start positions}})}, {Helper Functions: Build history, Interpolate, NDSolve}, { , Plot}, { , NDSolve}, { , change one car's behavior}})

Car

Class[Car,Import:>{firstCar,PositionAndVelocity,t}] := {
name=$.z.ClassNumber,
Static.previous=Null,
mass=0,speedup=0,slowdown=0,breaklimits=0,δ=0,
history = {},equation=True,

Car[m_,accel_,break_,lims_,delay_]:= With[{prev=previous,This=This},
mass=m;speedup=accel;slowdown=break;breaklimits=lims;δ=delay;previous = This;
equation:=With[{mass=mass,speedup=speedup,slowdown=slowdown,δ=δ},If[prev===Null,
name[t]==If[t≥0,firstCar[t]],
mass * name''[t]==If[t≥0,Block[{prevpos,prevvel,pos,vel,dist,td=t-δ},
{prevpos,prevvel}=PositionAndVelocity[td,prev.history];
{pos,vel}=PositionAndVelocity[td,This.history];
dist = prevpos-pos;
Which[
dist<1.75vel,-Clip[slowdown*((Max[vel-prevvel,1])/(Max[dist,0.1])),breaklimits],
dist>3.5vel,speedup,
True,0
]
]]
]]
]
};

PositionAndVelocity=Compile[{t,{history,_Real,2}},Block[{slot,frac},
slot=0;
While[t>history[[slot+1,1]],++slot];
frac=(t-history[[slot,1]])/(history[[slot+1,1]]-history[[slot,1]]);
history[[slot,{2,3}]]+frac*(history[[slot+1,{2,3}]]-history[[slot,{2,3}]])
]];

Add Cache

Class[DelayedCar,Car, Import :> {t,τ,AddPosition}] := {
Static.cached={},
DelayedCar[m_,accel_,break_,lims_,delay_]:=(
Super[m,accel,break,lims,delay];
AppendTo[cached,This]
),
Static.cache[t_, vals_] := If[Not[NumberQ[t]],
Map[(#1.history = {Flatten[{0,#2}]})&[Sequence @@ #]&, Transpose[{cached, vals}]],
Map[(
#1.history = AddPosition[t,#2,#1.history]
)&[Sequence @@ #]&,Transpose[{cached, vals}]]
],
Static.HoldFirst.reset[steps_,objs_,pos_,start_,first_]:=(
cache[Remove, Thread[{pos,D[first[t],t]/.t->0}]];
Do[cache[t,(#.name)[0]+D[(#.name)[τ],τ]t/.τ->0/.start&/@objs],{t,-3.,-.1,.1}];
steps=0
)
};

AddPosition=Compile[{t,newpos,{history,_Real,2}},Block[{pos,droppos,tp,tn},
droppos=0;
While[t-history[[droppos+1,1]]>3,++droppos];(*use only last 3 secs*)
pos=Drop[history,droppos];
If[t<0,
Insert[pos,{t,newpos,Part[history,-1,3]},-2],
pos=If[t==Part[history,-1,1],
ReplacePart[pos,{t,newpos,Part[history,-1,3]},-1],
Append[pos,{t,newpos,Part[history,-1,3]}]
];
tp=pos[[-2,1]]-pos[[-3,1]];
tn=pos[[-1,1]]-pos[[-2,1]];
pos[[-2,3]]=((pos[[-2,2]]-pos[[-3,2]])/tp+(pos[[-1,2]]-pos[[-2,2]])/tn)/2;
pos[[-1,3]]=pos[[-2,3]]+tn/tp*(pos[[-2,3]]-pos[[-3,3]]);
pos
]
]];

Create Identical Cars

DelayedCar.cached={};
cars = Table[New.DelayedCar[1000, 1500, 10000, {100,10000}, 0.5], {i, 5}];

Start Conditions

startPos=Range[Length[cars]-1,0,-1]*55;
firstCar=Function[First[startPos]+30*#-(75/π)*((Tanh[#/2])^5)*Sin[2π((#)/40)]];
funs=Through[(cars.name)[t]];
start=Join[Thread[funs==startPos],Thread[D[funs,t]==D[firstCar[t],t]]]/.t->0

Solve

solve:=Block[{steps,result},
DelayedCar.reset[steps,cars,startPos,Rule@@@start,firstCar];
result=Timing[First[With[{funs = funs},
 sol = First[NDSolve[Join[cars.equation, start], funs, {t, 0, 200},
MaxStepSize ->1/10,
PrecisionGoal->6,
Method -> {"EventLocator", "Event" :> (
++steps;
DelayedCar.cache[t, funs];
endTime = t;
minDist = Min[Subtract @@@ Partition[funs, 2, 1]];
If[minDist <= 0.1,EmitSound[Sound[SoundNote[#,0.01, "Gunshot"]&/@Range[60]] ]];
minDist > 0.1
)}
]]
]]];
Print[{Last[result],"..."}," in ",First[result]," seconds and ",steps," steps "];
]
solve

Plot

plot:=Plot[Evaluate[funs-(D[firstCar[t],t]/.t->0)*t /. sol], {t, 0,endTime-0.1}, Frame -> True, Axes -> False, PlotRange -> {{0,200},{-200,400}}, ImageSize->600,PlotStyle -> Table[Hue[i/Length[cars]], {i, Length[cars]}], PlotLabel->DisplayForm[StyleBox["Position Relative to Average Speed","Subsubsection"]]]
plot

Model Variations

Slow Reaction

cars[[2]].δ=1.2;
solve;
cars[[2]].δ=0.5;
plot

Motorbike

cars[[2]].mass=250;
solve;
cars[[2]].mass=1000;
plot

Nervous Driver

cars[[2]].speedup = 4000;
cars[[2]].breaklimits={6000,10000};
solve;
cars[[2]].speedup = 1500;
cars[[2]].breaklimits={100,10000};
plot

({{Analog Computer Circuits, one class template -> many objects}, { , incremental development -> bottom up design}})

Idea

({{Class, Object}, {out=fun(in), create elements}, {connection, connect elements}, {automatic naming, visually check layout}, {collect equations, collect equations}, {graphics, NDSolve}, {specific behavior,  }})

Basic Element

Class[BaseModelElement] := {
name := $.Class,
Final.var[slot_:1] := $.x.slot,
Final.call[slot_:1] := name[var[slot]],

sockets = {},
BaseModelElement[slots_Integer] := sockets = Table[Null, {slots}],
Final.plugIn[element_, slot_:1] := sockets = ReplacePart[sockets, element, slot],

Final.in[] := call/@Range[Length[sockets]],
Abstract.Virtual.Protected.out[] := $Failed,

Final.equations[] := Thread[in[]==sockets.out[]],
Static.Final.equations[elements__] := Module[{eqns},
eqns = Join@@(Flatten[{elements}].equations[]);
FixedPoint[Function[eq,
Map[D[#,Union[Extract[eq,Append[#,2]&/@Position[eq,_Integrate]]]]&,eq]
]/@#&,eqns] (*Integrate -> Derivative*)
]
};

Add Graphics

Class[ModelElement, BaseModelElement] := {
	inPos = {}, outPos = {}, back = False, rect = {},
	ModelElement[pos__, pt_, dir_] := (
		Super[Length[{pos}]-1];
		inPos = pt+#&/@If[back = dir===Backward, Reverse, Identity][{pos}];
		outPos = Last[inPos]; 
		inPos = Most[inPos];
		rect = {{GrayLevel[0.6], Rectangle[pt,pt+1]}, Text[name,pt+1/2], 
			{PointSize[0.015], Point/@inPos},{PointSize[0.03], Point[outPos]}}; 
	),
	graph[] := MapIndexed[Block[{fr = #.outPos, to = Extract[inPos, #2]}, 
		{rect,Line[{fr,Extract[If[#.back,{to,fr},{fr,to}],{{1,1},{2,2}}],to}]}
	]&, sockets]
};

Specific Elements

Integrator

Class[Integrator,ModelElement] := {
Integrator[x_, y_, dir_:False] := Super[{0,1/2},{1,1/2}, {x, y}, dir], 
Override.Virtual.out[] := Integrate[call[1], var[1]]
};

Gain

Class[Gain,ModelElement] := {
Private.factor = 1,
Gain[x_, y_, dir_:False, factor_] := 
(Super[{0,1/2},{1,1/2}, {x, y}, dir]; This.factor = factor), 
Override.Virtual.out[] := factor*call[1]
};

Feedback

Class[Adder,ModelElement] := {
	Adder[x_, y_, dir_:False] := Super[{0,1/2},{1/2,1},{1,1/2}, {x, y}, dir], 
	Override.Virtual.out[] := Sum[call[i], {i,Length[sockets]}]
};

Objects

amp1 = New.Gain[2,3,Backward,-a];
int1 = New.Integrator[1,0];
int2 = New.Integrator[3,0];
objs={int1, int2, amp1};

Connections

int1.plugIn[amp1];
int2.plugIn[int1];
amp1.plugIn[int2];

Show[Graphics[objs.graph[]], AspectRatio -> 1/2, TextStyle -> "Subsubtitle"]

Equations

funs=Flatten[objs.sockets.in[]]
vars=Union[Sequence@@@funs]
eqns=ModelElement.equations[objs]

start={(int2.name)[0]==0,int2.name'[0]==1}

Solutions

DSolve[Join[eqns,start],funs,vars]

Add Damping

add1 = New.Adder[-1,0];
amp2 = New.Gain[0,1.5,Backward,-d];
objs={int1, amp1, int2, amp2, add1};

amp2.plugIn[int1];
add1.plugIn[amp1];
add1.plugIn[amp2, 2];
int1.plugIn[add1];

Show[Graphics[objs.graph[]], AspectRatio->1/2, TextStyle->"Subsubtitle",ImageSize->600]

funs=DeleteCases[Flatten[objs.sockets.in[]],(add1.name)[_]]
vars=Union[Sequence@@@funs]
eqns=ModelElement.equations[objs];
Apply[Rule,Select[eqns,Not[FreeQ[#,(add1.name)[_]==_]]&],{1}];
eqns=DeleteCases[eqns//.%,True]

FullSimplify[DSolve[Join[eqns,start],funs,vars]]/.(d^2)->Δ+4a//Timing

({{Further Potentials,  }, { ,  }})

({{Geographic Information System, masking objects,layerclasses}, {Ray Tracing, lenses,rays}, {Multi-Body Simulation, nonuniform systems}, {Supply Chains, queues}, {Portfolio Optimization, individual greeks}, {Reaction Kinetics, molecule classes}})

({{Objectica, }})

({{({{Industry, 1000$+tax, single user, Ringstraße 100}, {Education, 1000 $+tax, staff (max 30)+students, 52078 Aachen}, {Professorship, 666 $+tax, staff (max 30), Germany}, {University, 333 $+tax, single user, Tel +49/241/566 08 28}, {Student, 33 $+tax, single user, sales@symbols-and-numbers.com}})}, {→ Special Conference Offer:-20% ←}})

({{abstraction layer, stepwise refinement -> iterative model development}, {OO libraries, black box models -> nonexpert can build models}, {one class template -> many objects, easy to change -> refinement driven by simulation results}, {encapsulation of functionality, incremental development -> bottom up design}})