A Smallest Example for Multithreading in Python with Event and Lock
Thanks to Jacob Zhong’s help. This is quite useful!
1 

Transportation Research Board 2019
Topic:
 What proportion of the vehicles with known trajectories is sufficient to estimate the complete information?
Journal Club 02 Introduction to Convex Optimization
The slides are from the Journal club 02.
Questions remain to be solved:
 What is the difference between p1 norm and the p2 norm in the linear regression weight decay?
 Why the dual Lagrangian problem always convex?
 Need to know more about the VCdimension.
Realtime Signal Control System using Trajectories
Remain to be done!
2018 Christmas Trip in Los Angeles, Las Vegas, Grand Canyon and Antelope Canyon
2018 Christmas trip with Laixi, Ruxing, Zihan from 23th to 27th Decemember, 2018.
The complete version of the video:
If you can not get access to YouTube, you can find the video in Bilibili [here]:
Day 1: Disneyland
Day 2: From LA to Las Vegas
Day 3: Grand Canyon
Day 4: Antelope Canyon
CEE 553 Infrastructure System Optimization Finalterm Exam Review
This post is the continued part of [link].
Constrained Optimization
Interpretation of Lagrangian Multipliers
Consider the family of problems:
$$\min f(x)  h(x)=u.$$
For each $u$, $x^*(u)$ and $\lambda ^*$ are the corresponding regular local minimum and Lagrangian multiplier respectively. Denote $p(u)=f(x^*(u))$, then:
$$\nabla p(u) =  \lambda ^* (u)$$
where the $\lambda ^*$ can be interpreted as shadow price.
Here gives a brief proof:
Assume the original constrained optimization problem is: $$\min f(x), s.t.\quad h(x)=u,u=Const$$ then according to the KKT condition: $$\nabla f(x(u))+\lambda \nabla h(x(u))=0$$ $$ \frac{dx(u)^T}{du}(\nabla f(x(u))+\lambda \nabla h(x(u)))=0$$ where for the $h(x(u))$: $$h(x(u))=u\rightarrow \frac{dh(x(u))}{du}=\frac{dx(u)^T}{du}\cdot \nabla h(x(u)) =1$$ then finally: $$\frac{f(x^*(u))}{du}=\frac{dx(u)^T}{du}\cdot \nabla f(x(u))=\lambda$$
Solving Constrained Optimization
Onedimension Condition
For the onedimension constrained optimization problem,
$$\min_{x\in [a,b]}z(x)$$
where $z(x)$ is convex and continuous. Here are two methods: bisection method and golden section method:
Bisection Method
 Let $[a^n, b^n]$ be the current interval
 $\tilde{x}^n=\frac{a^n+b^n}{2}$
 If $\frac{dz(\tilde{x}^n)}{dx}<0$, discard $[a^n,\tilde{x}^n]$ and $\tilde{x}^n,b^n \rightarrow a^{n+1},b^{n+1}$; otherwise discard $[\tilde{x}^n,b^n]$ and $a^n,\tilde{x}^n \rightarrow a^{n+1},b^{n+1}$;
 If $b^{N}a^N<\epsilon$, return. Otherwise go to step 2.
Golden section method is similar to the bisection method while the main difference is that you do not have to calculate the derivative of the lost function and for each iteration you need only calculate one extra value.
Golden Section Method
 Let $[a^n,b^n]$ be the current interval
 Set $x_L^n=0.618a^n+(10.618)b^n$ and $x_R^n=(10.618)a^n+0.618b^n$
 If $z(x_L^n)<z(x_R^n)$, discard $[x_R^n,b^n]$ and $a^n,x_R^n\rightarrow a^{n+1}, b^{n+1}$; otherwise discard $[a^n,x_L^n]$ and $x_L^n,b^n \rightarrow a^{n+1}, b^{n+1}$;
 If stopping criterion met, return; otherwise go step 2.
Normal Situation
By applying the thought of the iterative descent method to the constrained optimization problem, the feasible descent direction is:
$$\nabla f(x^k)^T(yx^k)<0$$
where $y\in X$ and $X$ is a convex region.
The step size can be chosen by searching for a minimum between $y$ and $x^k$.
Conditional Gradient Method (FrankWolfe Method)
For a linear constrained problem:
$$\min_{x\in R^n}z(x)\quad s.t.\sum_i h_{ij} x_i \ge b_j\forall j\in J$$
 Start with a feasible solution $x^0$, $n=1$;
 Direction finding. Find $y^n$ that solves the linear program:
 $$\min z_L^n(y)=\sum_i (\frac{\partial z(x^n)}{\partial x_i})y_i, s.t. \sum_i h_{ij}y_i\ge b_j \forall j\in J$$
 Step size determination. Find $\alpha_n$ that solves:
 $$\min_{0\le \alpha \le 1}z[x^n+\alpha (y^nx^n)]$$
 Move. Set $x^{n+1}=x^n+\alpha_n(y_nx_n)$.
 Convergence test. If $z(x^n)z(x^{n+1})\le \kappa$, stop. Otherwise, let $n=n+1$ and go to step 1.
Another method called Gradient Projection Method is first find the gradient descent direction and then project the direction to the feasible plate.
Linear Programming
Standard Form of LP
The standard form of the linear programming is:
$$\min c^T x, s.t.Ax=b,x\ge 0$$
so the standard linear programming problem is convex program because both the function and the feasible region are convex. Then the KKT conditions are both necessary and sufficient.
There are only three possible results for a linear programming problem:
 Has an optimal solution (may be multiple global optimum solutions)
 Has no optimal solution because the problem is not bounded
 No feasible region.
Primal and Dual Problems
The primal problem:
$$\min c^T x ,s.t.Ax=0,x\ge 0$$
and the dual problem:
$$\max p^Tb,s.t.A^Tp\le c$$
They are called primal problem and the dual problem because they share the same KKT conditions:
$$c^Tx=p^Tb,Ax=b,A^Tp\le c,x\ge 0$$
If both $x$ and $p$ are feasible, then we have: $c^Tx\ge p^Tb$. At optimality, $c^Tx^*=p^{*T}b$. There are two special conditions: if the primal is unbounded, then the dual is infeasible; if the primal is not feasible, then the dual will be unbounded.
Simplex Method
General procedure of simplex method
 Choose a starting feasible basis, B.
 Compute $B^{1}$, then the basic feasible solution is $x_B=B^{1}$ and $x_N=0$ and the objective value is $c^Tx=c^T_BB^{1}b$.
 Calculate the reducted cost for each nonbasic variable, $\bar{c}_j=c_jc^T_BB^{1}A_j$. If all reducted costs are greater or equal than zero, stop and the current basic feasible solution is optimal. Otherwise, select a nonbasic variable, $x_j$, with a negative reduced cost and to to step 3.
Calculate $d_B=B^{1}A_j$. If $d_B\ge 0$, stop and the problem is unbounded. Otherwise, let $x_{B(l)}$ be the variable which yields the maximum feasible step size,
$$\frac{x_{B(l)}}{d_{B(l)}}=\min (\frac{x_{B(i)}}{d_{B(i)}})$$
Replace the column$A_{B(l)}$ with column $A_j$ to obtain a new basis, $B$. Go to step 1.
Here is a brief interpretation of simplex method:
First, we select $B$ as the basis matrix and the corresponding variables are called basic variables, which means that these variables do not have to be zero (they may be zero, they just do not have to be). Then the leftover variables are called nonbasic variables which are all equal to zero. When we selct the matrix $B$ and if $B$ is invertible, we can get a series value of the basic varialbes. Besides, all the value of the basic variables have to be greater than zero otherwise they will not satisfied the constraints. So usually finding the starting feasible basis is not that easy and later will discuss about how to find the feasible basis.
Assume that we have found the basis matrix $B$ and also the feasible basic variables. Then we can move from this point to another joint vertex by choosing a entering variable and a leaving variable. Then along the direction of each entering variable, the reduced cost can be calculated. If all the feasible directions towards a candidate entering variables are not the reduced direction, then this point is the optimal point. Otherwise, we can always find a reduced direction, which means that we can determine an entering variable. After entering variable we can move towards this entering variable and see which basic variables become zero first. If no variable is reduced along the direction, then this problem will be unbounded, otherwise, we can find a leaving variable. Then we have finished an iteration to the next point.
Finding the start feasible basis
Sometimes it is hard to find the start feasible basis. One possible way to find a feasible basis is to add extra variables to be the basic variables and then set their weights as large numbers.
Sensitivity Analysis
 If the cost vector $c$ changes, then the current basis will still be the optimal if the reduced costs are greater zero for all nonbasic variables, otherwise continue the iteration of simplex method. Things are similar if the matrix $A$ changes.
If you have a better understanding towards the simplex method, it is not difficult to analyze the possible result if we change some of the parameters.
Code of Numerical Methods
I have coded most of the numerical methods included in this course and the code has been pushed to GitHub, see link.
2018 Michigan Traffic Lab ThanksGiving Day Dinner at Dr. Feng's Home
感恩节晚上我们整个实验室在冯毅恒老师家里组织了感恩节晚餐。
大家聚在一起一起吃火锅，玩游戏，可以说是相当热闹了。
Hidden Markov Map Matching  Method and Experiment
Reference
This essay is based on the following paper:
 Newson, Paul, and John Krumm. “Hidden Markov map matching through noise and sparseness.” Proceedings of the 17th ACM SIGSPATIAL international conference on advances in geographic information systems. ACM, 2009. [pdf]
Introduction
Matching GPS points to roads is the foundation of all kinds of traffic measurement using trajectories data as the input. This essay will give a brief introduction to the trajectories map matching method and I also write a project based on Python to use this method to match the data.
Here is a description of the map matching problem, the grey lines in the figure are the road of Ann Arbor while the light blue lines are trajectories data (GPS points). However, we cannot determine which road the GPS point belonged to without map matching and this is what map matching can do.
I have put the project in my GitLab account, this project is open for everyone in Michigan Traffic Lab, feel free to request the right to edit/view this project if you are a member of Michigan Traffic Lab or a student in U of M. I will be glad to share with you if you are interested. [Project Link Here]
Method
Model Establish
The key problem in map matching is the tradeoff between the roads suggested by the location data and the feasibility of the path. One simple method to determine the link the one GPS point belongs to is the find the nearest link among all the candidate links. However, this method ignore the relation between the GPS points and the feasibility of the path.
The following figures give the method in the paper to do the map matching using a Hidden Markov Model. (Figures are also from the reference)
There are two parts of the model, the first is the measurement probabilities:
$$p(z_tr_i)=\frac{1}{\sqrt{2\pi}\sigma_z}e^{0.5(\frac{z_t x_{t,i}_{greatCircle}}{\sigma_z})^2}$$
Then is the transition probabilities:
$$p(d_t)=\frac{1}{\beta}e^{\frac{d_t}{\beta}}$$
where the $d_t$ equals to:
$$d_t=z_tz_{t+1}_{greatCircle}x_{t,i^*}x_{t+1,j^*}_{route}$$
Then this problem can be converted to a dynamic programming problem (DP).
Dynamic Programming
Some references here:
Experiment
The input of the algorithm is the map information and the trajectories data shown as below.
To match the trajectory to the map, we first determine the block zone of the trajectory (red block in the figure) and then remove the road that are away from the block zone. Then we can calculate the distance between the trajectory point and the road. The a threshold was chosen to roughly select the candidate roads. As shown in the figure, the blue points are the trajectory points and the green points are the candidiate projected points. The red points are the mapmatching results using dynamic programming.
The following figure gives the result of the dynamic programming. Axis x is index of time stamp (from $t=0$) and axis y is the index of the candidate road. The black lines are the whole network and the red line shows the global optimal solution of the problem (the chosen road).
Here are some mapmatching results using the trajectory data with sampled period $T=60s$.